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Preface 


This  study  was  designed  to  research  and  develop  the  methodology  necessary 
to  implement  a  general  computer  program  that  determines  an  approximate  solu¬ 
tion  to  a  system  of  Fredholm  integral  equations.  The  computer  program  was 
designed  to  be  easily  converted  into  an  excellent  research  and  analysis  tool  for  the 
Defense  Nuclear  Agency’s  underground  nuclear  effects  simulation  testing  program. 

Due  to  the  symbols  used  during  the  study,  two  lists  of  symbols  are  provided. 
The  first  List  of  Symbols  contains  the  symbols  used  to  develop  the  methodology 
for  this  study.  The  second  list  is  located  in  Tab  1  of  Appendix  C  and  contains  all 
the  variables  used  in  the  development  of  the  computer  program. 

In  conducting  this  study.  I  received  guidance  and  strength  from  a  number  of 
people.  First.  I  would  like  to  thank  my  thesis  advisor  LCDR.  Kirk  A.  Mathews 
for  the  assistance  and  knowledge  he  provided  during  my  study.  I  would  also  like 
to  thank  my  thesis  committee;  Dr.  George  John,  Maj.  Jim  A.  Lupo,  and  Lt.  Col. 
Albert  Alexander  (from  the  Defense  Nuclear  Agency)  for  their  assistance  in  com¬ 
piling  and  completing  this  final  copy  of  my  thesis.  Also  a  special  thanks  to 
LCDR.  Mathews  and  Lt.  Col.  Alexander  for  their  assistance  in  arranging  my 
assignment  to  Field  Command  Defense  Nuclear  Agency  as  radiation  diagnostician 
based  on  the  knowledge  gained  during  this  study.  Lastly.  I  would  like  to  thank 
my  wife,  Karen,  for  her  needed  assistance  in  typing  as  well  as  moral  support  dur¬ 
ing  this  study. 
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Abstract 


The  purpose  of  this  study  was  to  develop  the  methodology  for  and  to  imple¬ 
ment  a  computer  program  to  approximate  a  solution  to  a  system  of  Fredholm 
integral  equations.  The  system  of  equations  used  in  this  study  is  representative  of 
the  equations  formed  during  the  detection  of  pulsed  radiation  using  a  series  of 
detectors  with  asymmetric  response  functions.  Though  general  in  nature  and 
applicable  to  all  systems  of  Fredholm  integral  equations,  the  equations  studied  are 
of  importance  to  the  Defense  Nuclear  Agency  with  regard  to  the  measurement  of 
radiation  spectra  during  underground  nuclear  effects  simulation  testing. 


The  deconvolution  or  solution  technique  consisted  of  representing  the 
unfolded  spectrum  as  a  weighted  sum  of  basis  functions.  This  unfolded  spectrum, 
the  actual  spectrum,  and  a  predicted  spectrum  were  then  used  to  form  a  x2  test 
statistic.  By  adjusting  the  parameters  in  the  basis  functions  and  their  weights.  \2 
was  minimized  and  the  unfolded  spectrum  was  corrected  to  approximate  the 
actual  spectrum. 


The  methodology  for  this  deconvolution  technique  was  then  converted  into  a 

general  computer  program.  The  validation  cases  conducted  on  the  two  types  of 

spectra  confirmed  the  reliability  of  the  methodology  and  the  computer  program^ 

Additionally,  an  initial  study  with  simulated  measurement  error  added  to  the 

measured-to-predicted  ratios  showed  that  the  actual  spectrum  could  not  be 

returned  exactly.  The  second  study  approximated  the  actual  spectrum  with  an 

S 

unfolded  spectrum  using  a  second  set  of  basis  functions.  An  acceptable  approxi¬ 
mation  was  conducted;  however,  certain  artifacts  were  discovered  in  the  unfolded 
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The  validation  cases  and  preliminary  test  cases  conducted  prove  that  the 
computer  program  based  on  the  methodology  presented  in  this  study  is  a  viable 
means  of  approximating  an  actual  radiation  spectrum.  Using  this  study  and  com¬ 
puter  program  as  a  starting  point,  the  study  of  new  basis  functions  and  the  effect 
of  how  well  the  actual  spectrum  can  be  approximated  based  on  the  number  of 
detectors  available  to  determine  the  spectrum  is  recommended. 
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AN  APPROXIMATION  TECHNIQUE  FOR  SOLVING 
A  SYSTEM  OF  FREDHOLM  INTEGRAL  EQUATIONS 
FOR  ASYMMETRIC  DETECTOR  RESPONSE  FUNCTIONS 


It  Introduction 

Background 

Approximate  solutions  to  systems  of  Fredholm  integral  equations  are  needed 
by  the  Air  Force  in  many  fields  of  study,  ranging  from  acoustics  to  optics.  Such 
solutions  also  occur  in  disciplines  as  varied  as  geology  and  astronomy.  An 
instance  of  particular  importance  to  the  Air  Force  arises  in  the  measurement  of 
radiation  spectra  emitted  by  nuclear  devices  detonated  underground  by  the 
Defense  Nuclear  Agency  for  simulation  of  nuclear  weapons  effects.  This  study 
develop^  an  approximate  solution  for  a  set  of  Fredholm  integral  equations  of  type 
1  of  the  following  form: 

OO 

Y,m  -  JS.(E)R,(E)dE  (1) 

o 

where 

Y,m  *=  the  measured  signal  of  the  ith  detector 
Sa(E)  ■*  the  actual  spectrum  as  a  function  of  energy 
R,  (E )  =  the  calibrated  response  function  of  the  ith 
detector  as  a  function  of  energy 

In  the  detection  of  pulsed  radiation,  a  set  of  detectors  covering  various  energy 
ranges  is  used  to  obtain  a  set  of  measured  signals.  Each  of  these  measured  signals 
can  be  represented  by  an  equation  of  the  form  of  Equation  (1).  Since  the  response 
functions  of  these  detectors  are  not  rectangular  in  shape  with  neglibie  width  and 
since  a  finite  number  of  detectors  must  be  used,  only  a  limited  resolution  of  the 
actual  spectrum  can  be  achieved.  In  order  to  achieve  a  reasonable  resolution  a 
filter-fluorescer  detection  system  is  utilized.  The  detectors  and  their  response 
functions  are  discussed  in  detail  in  Section  IT. 


Also,  the  experimenter  does  not  know  the  actual  spectrum  or  the  exact 
response  functions  of  the  detectors.  However,  the  experimenter  can  calibrate  the 
detectors  used  and  produce  calibrated  response  functions  for  the  detectors.  The 
experimenter  can  also  predict  the  shape  of  the  actual  spectrum  based  on  source 
design  and  previous  measurements.  Thus,  at  the  end  of  the  experiment,  the 
experimenter  has  the  measured  signals,  a  predicted  spectrum,  and  the  calibrated 
response  functions  of  the  detectors  to  use  in  approximating  the  actual  spectrum. 
The  procedure  used  to  conduct  this  approximation  is  called  deconvolution  or 
unfolding. 

Problem 

The  measured  signals  discussed  above  are  the  best  information  the  experi¬ 
menter  has  concerning  the  ideal  signal  (the  error  free  signals  from  the  actual  spec¬ 
trum.)  These  measured  signals  contain  recording  error  and  transmission  error  in 
addition  to  the  error  due  to  the  response  functions.  These  errors  are  discussed  in 
detail  in  Section  II.  The  main  problem  addressed  in  this  study  was  the  develop¬ 
ment  of  the  methodology  necessary  to  conduct  a  general  deconvolution  of  the 
actual  spectrum  from  the  set  of  calibrated  response  functions  and  the  set  of  meas¬ 
ured  signals  noted  in  Equation  (1).  In  other  words,  the  known  measured  signals, 
calibrated  response  functions,  and  the  predicted  spectrum  are  used  to  determine 
an  unfolded  spectrum.  When  this  unfolded  spectrum  is  folded  with  the  calibrated 
response  functions,  an  approximation  to  the  measured  signals  are  returned.  This 
unfolded  spectrum  is  then  used  as  an  approximation  to  the  actual  spectrum. 
Secondly,  a  computer  program  was  developed  to  conduct  this  deconvolution  using 
various  basis  functions  to  construct  the  unfolded  spectrum. 


Since  unclassified  experimental  data  was  unavailable,  this  study  considered 
spectra  that  were  constructed  of  either  normalized  planckian  black  body  distribu¬ 
tions,  called  planckians,  or  of  cubic  splines.  During  the  study,  a  trial  spectrum 
was  constructed  and  used  as  the  actual  spectrum.  This  spectrum  was  then  used 
to  calculate  the  ideal  signals  or  measured  signals  if  simulated  measurement  errors 
are  introduced.  First,  a  study  was  conducted  to  determine  if  the  basis  functions 
could  be  used  to  unfold  an  actual  spectrum  produced  from  the  same  type  of  basis 
functions  with  no  error  introduced.  Then  a  study  was  conducted  to  demonstrate 
the  effects  of  simulated  measurement  error  in  the  measured-to-predicted  ratios  on 
the  unfolded  spectrum.  This  simulated  error  represented  the  errors  between  the 
ideal  signal  and  the  measured  signal.  Finally,  a  study  was  conducted  to  deter¬ 
mine  if  the  basis  functions  could  be  used  to  unfold  an  actual  spectrum  constructed 
from  another  set  of  basis  functions.  The  energy  range  considered  in  this  study 
was  from  E ,°  to  128£,°,  where  E?  (the  k  edge  of  the  fluorescer  for  the  first  detec¬ 
tor)  is  used  as  a  convenient  arbitrary  energy  unit. 

Assumptions 

In  the  development  of  the  methodology,  a  basic  assumption  was  that  the 
nonuniqueness  of  the  actual  spectrum  and  the  errors  developed  in  the  mathemat¬ 
ics  of  an  analytic  solution  required  the  approximation  technique  to  be  numerical. 
Also  in  order  to  define  or  bound  the  study  and  to  allow  a  test  method  to  be 
developed,  the  following  assumptions  are  made: 

1.  Twenty  detectors  are  available  to  conduct  the  experiment.  Thirteen 
detectors  have  a  closed  response  function  and  seven  detectors  have  an 
open  response  function,  as  described  in  Section  III. 

2.  A  resolving  power  (E  center/  A  E)  of  about  1.5  is  desired. 


3.  The  exact  response  functions  are  equal  to  the  calibrated  response  func¬ 
tions,  as  described  in  Section  II.  (This  does  not  restrict  the  applicability 
of  this  analysis,  since  calibration  errors  are  indistinguishable  from  other 
measurement  and  recording  errors.) 

General  Approach 

The  approach  used  in  this  study  was  to  define  an  unfolded  spectrum  using  a 
series  of  basis  functions: 

S.iE)=£<‘,r,{E.P)  (2) 

y-i 

where 

SU(E)  =  the  unfolded  spectrum 

a;  =  the  coefficient  for  the  jth  basis  function 
F j[E ,P)  =*  the  jth  basis  function 

nj  =  the  number  of  basis  functions 
E  =  energy 

P  =  the  parameter(s)  of  F; 

The  next  step  was  to  compute  a  v2  test  statistic  using  the  ratios  of  the  ideal 
signal  to  the  predicted  signal  and  the  ratios  of  the  unfolded  signal  to  the  predicted 
signal,  together  with  the  measurement  uncertainties  of  these  ratios.  These  ratios 
are  defined  in  detail  in  Section  III.  This  \2  test  statistic  was  then  minimized  using 
the  Fletcher-Powell  technique  to  vary  the  parameters  and  weighting  factors  or 
coefficients  in  the  basis  functions  of  the  unfolded  spectrum.  The  term  parameters 
will  refer  to  the  parameters  of  the  function  as  well  as  the  coefficients  for  the  func¬ 
tion  in  the  remainder  of  this  study. 

Sequence  Of  Presentation 

Section  II  of  this  study  presents  a  detailed  analysis  of  the  problem  including 
other  deconvolution  methods  currently  being  used  and  possible  errors  to  be  con- 


sidered.  Section  III  presents  the  theory  used  to  develop  the  computer  program 
while  Section  IV  presents  the  development  of  the  computer  program.  Section  V 
presents  the  results  and  a  discussion  of  these  results  for  the  validation  cases  and 
the  test  cases  conducted  using  the  computer  program.  The  test  cases  investigated 
the  strengths  and  weaknesses  of  the  two  types  of  basis  functions  and  the  effect 
simulated  measurement  error  in  the  measured-to-predicted  ratios  had  on  the 
unfolding  technique  using  these  basis  functions.  Section  VI  then  summarizes  the 
study  and  presents  my  conclusions.  Recommendations  for  future  studies  are 
presented  in  Section  MI. 


II.  Detailed  Problem  Analysis 


Introduction 

In  Section  I,  the  main  problem  addressed  was  stated  to  be  the  determination 
of  the  methodology  and  the  implementation  of  a  computer  code  to  approximate 
the  actual  spectrum  in  Equation  (1),  the  Fredholm  integral  equation.  The  pur¬ 
pose  of  this  section  is  to  describe  the  types  of  errors  found  in  unfolding  tech¬ 
niques,  to  explain  how  these  errors  were  handled,  and  to  present  the  two  most 
popular  unfold  techniques. 

Errors  Introduced  During  Unfolding 

During  the  unfolding  process,  four  main  errors  contribute  to  the  uncertainty 
in  the  unfolded  spectrum.  The  6rst  is  simply  the  measurement  error  introduced 
in  the  measured  signals.  If  the  same  test  was  conducted  a  number  of  times,  a 
mean  value  for  the  measured  signals  could  be  established  along  with  a  standard 
deviation.  This  measurement  error  includes  the  errors  introduced  form  the 
transmission  of  the  ideal  signals,  the  recording  of  the  ideal  signals,  and  the  read¬ 
ing  of  the  ideal  signals.  However,  in  underground  testing  and  the  detection  of 
cosmic  radiation  by  satellites,  the  experimenter  only  makes  one  test.  Thus,  the 
experimenter  must  be  aware  of  and  account  for  the  possible  error  introduced  by 
the  statistical  nature  of  the  measured  signals. 

The  second  error  is  found  in  the  response  functions  of  the  detectors.  No 
matter  how  carefully  an  experimenter  calibrates  the  detectors,  one  will  not  be  abl 
to  determine  the  exact  response  functions  of  the  detectors.  This  difference 
between  the  exact  response  functions  and  the  calibrated  response  functions  is  the 
second  type  of  error  introduced  in  the  unfolding  procedure. 

The  third  error  is  found  in  the  mathematics  of  the  unfolding  process.  The 
error  is  introduced  by  converting  the  Fredholm  integral  into  a  summation  over 


very  small  energy  intervals.  However,  this  error  can  be  reduced  so  as  to  be  negli¬ 
gible  compared  to  the  other  errors  by  selecting  the  proper  energy  intervals  or  bin 
widths  to  evaluate  the  integral. 

The  final  error  to  be  discussed  is  that  caused  by  the  nonuniqueness  of  the 
solution  or  approximation  to  the  actual  spectrum.  Because  the  approximation  of 
the  actual  spectrum  requires  an  infinite  amount  of  detail  to  be  resolved  from  a 
finite  amount  of  information,  the  solution  is  non-unique  and  the  problem  is  ill- 
posed.  The  limited  number  of  detectors  and  their  response  functions  ensure  that 
the  experimenter  can  not  resolve  all  the  detail  of  the  spectrum.  Thus,  the 
unfolded  signals  of  a  number  of  spectra  could  be  identical  even  though  the  spectra 
are  different.  If  a  detector  with  a  response  function  which  is  shaped  like  a  rectan¬ 
gle  with  neglible  width  was  available  at  every  energy  level,  the  actual  spectrum 
could  be  determined  exactly.  Since  this  is  not  possible  or  feasible,  the  experi¬ 
menter  must  be  reminded  of  the  nonuniqueness  of  the  approximation  to  actual 
spectrum  or  the  unfolded  spectrum. 

Errors  in  Actual  Practice  Versus  Simulated  Errors  in  This  Study 

Of  the  four  errors  introduced  above,  two  are  handled  differently  in  this  study 
than  they  are  in  actual  practice.  In  actual  practice,  the  measurement  error  and 
the  error  introduced  by  not  knowing  the  exact  response  functions  are  present  in 
the  measured  signals.  These  errors  are  then  carried  over  to  the  measured-tcv 
predicted  ratios.  However  in  order  to  account  for  these  errors,  the  experimenter 
will  approximate  the  standard  deviation  of  these  measured-lo-predicted  ratios,  rr_ , 
and  use  this  as  a  weighting  factor  for  the  test  statistics  as  discussed  in  the  next 
section.  This  <7,  is  based  on  the  knowledge  gained  from  past  experiments  and  the 
calibrations  conducted  on  the  detectors. 

In  this  study,  it  was  assumed  that  the  exact  response  functions  were  equal  to 
the  calibrated  response  functions  and  since  the  signals  were  not  measured,  recorde. 


or  transmitted,  the  measurement  errors  normally  found  in  the  measured  signals 
were  not  present.  Thus,  when  an  actual  spectrum  was  folded  with  the  calibrated 
response  functions  to  determine  a  set  of  measured  signals,  these  two  types  of 
error,  the  measurement  error  and  the  calibration  error,  were  not  present  and  the 
measured  signal  was  equal  to  the  ideal  signal.  In  order  to  introduce  these  errors 
in  this  study,  an  option  to  apply  a  normal  distribution  to  the  measured-to- 
predicted  ratios  was  included.  This  is  discussed  in  detail  in  Section  ID. 

The  third  error,  the  error  due  to  evaluating  the  Fredholm  integrals,  is  present 
in  both  the  actual  case  and  this  study.  As  discussed  in  Section  III,  the  error 
becomes  negligible  in  both  cases  by  selecting  a  small  bin  width.  Finally  the  fourth 
error  is  also  a  part  of  both  cases  and  was  demonstrated  by  studying  the  unfolding 
of  an  actual  spectrum  using  the  two  different  sets  of  basis  functions. 

The  mathematical  propagation  of  the  four  errors  discussed  above  was  not 
considered  in  this  study.  However,  the  propagation  of  the  errors  was  demon¬ 
strated  by  the  study  conducted  using  the  simulated  measurement  error.  Also,  the 
calibration  error  may  be  considered  as  part  of  the  simulated  measurement  error. 

Current  Unfolding  Techniques 

A  number  of  computer  programs  have  been  developed  to  approximate  the 
actual  spectrum  implicitly  given  in  Equation  (1).  However,  two  techniques,  the 
iterative  technique  and  the  cubic  spline  technique,  are  the  most  common  (6).  This 
section  will  describe  the  two  techniques  and  some  of  their  weaknesses. 

Iterative  Technique  (6).  The  widely  used  iterative  technique  starts  by 

folding  a  trial  spectrum  with  the  calibrated  response  functions  of  the  detectors 
and  comparing  these  measured  signals  to  the  measured  signals  from  the  experi¬ 
ment.  The  trial  spectrum  is  then  modified  and  smoothed  to  ensure  non-negativity 
and  continuity.  The  procedure  is  then  repeated  until  the  two  sets  of  measured 


'  VJf  V”  -  *>  ’>  ”  >'-•  V'J^>  ->  ■>  \>  ■>  V  V  V  V  ■>  W  •.“  V  V  W  W  -  J<  ■/  V  ’/*/  V  W  V ~S 


signals  converge.  By  allowing  the  measured  signals  to  converge,  various  artifacts 
such  as  discontinuities  and  spikes  are  usually  introduced  in  the  trial  spectrum. 
Thus,  the  procedure  is  generally  continued  only  until  the  computed  signals  for  the 
trial  spectrum  are  brought  into  an  acceptable  agreement  with  the  measured  sig¬ 
nals,  but  stopped  before  unacceptable  artifacts  develop  in  the  trial  spectrum. 

The  iterative  technique  has  three  main  limitations.  First,  because  the  itera¬ 
tion  is  not  allowed  to  converge,  the  shape  of  the  unfolded  spectrum  depends  on 
the  initial  guess  at  the  trial  spectrum.  Secondly,  the  technique  provides  little 
information  to  allow  for  error  analysis.  Finally,  the  decision  to  stop  the  iteration 
before  unacceptable  artifacts  develop  forces  the  technique  to  be  user  dependent. 

Cubic  Spline  Technique  (6).  The  cubic  spline  technique  used  here  is  based 
on  the  method  developed  at  the  Lockheed  Palo  Alto  Research  Laboratory.  The 
technique  consists  of  building  an  unfolded  spectrum  from  a  set  of  cubic  splines 
with  each  cubic  spline  in  the  set  being  multiplied  by  an  expansion  coefficient.  The 
cubic  splines  used  are  piecewise  I.agrangian  interpolating  splines  and  are  not  the 
common  B-splines.  In  the  cubic  spline  technique  the  knot  locations  are  adjusted 
and  the  number  of  knots  used  is  varied  to  obtain  the  best  unfolded  spectrum  or 
approximation  to  the  actual  spectrum.  The  cubic  spline  technique  is  an  excellent 
approximating  technique:  however,  the  main  limitation  of  the  cubic  ■spline  tech¬ 
nique  is  t  lie  dependence  on  the  knot  locations  or  points  used  to  develop  the  cubic 
'p!in**s.  Thi-  limitation  is  reduced  bv  selecting  the  best  locate  ns  f<  >r  the  knots 
before  developing  the  cubic  splines. 


III.  Theoretical  Development 

Introduction 

In  the  last  section,  the  overall  problem  and  the  two  most  common  techniques 
to  solve  this  problem  were  presented.  This  section  presents  the  theory  used  to 
develop  the  methodology  for  the  computer  program  in  this  study,  including  the 
selection  of  the  response  functions;  the  definition  of  the  various  spectra  and  sig¬ 
nals;  the  formulation  of  a  X'2  test  statistic;  the  minimization  of  this  x*  test  statis¬ 
tic;  the  calculation  of  a  normal  distribution  for  simulating  measurement  error;  and 
the  inclusion  of  a  flux  non-negativity  constraint. 

Response  Functions 

As  stated  in  Section  I,  this  study  assumes  the  exact  response  functions  are 
equal  to  the  calibrated  response  functions  and  that  twenty  detectors  will  be  used 
in  the  experiment.  Two  types  of  detector  systems  will  be  used  in  this  study.  The 
first  detection  system  is  a  filter-fluorescer  detection  system  and  produces  a  closed 
response  function.  This  closed  response  function  consists  of  a  section  that  is  simi¬ 
lar  to  a  narrow  rectangle  and  is  an  approximation  to  an  ideal  response  function. 
This  section  is  referred  to  as  the  '*inband"  response  and  has  a  response  between 
the  k-edge  of  the  fluorescer  and  the  k-edge  of  the  filter.  However,  a  second  section 
or  tail  section  is  also  present  in  the  closed  response  function.  The  tail  section  is 
formed  by  the  response  of  the  detector  to  energies  above  the  k-edge  of  the  filter. 
The  closed  response  function  is  depicted  in  Figure  1. 

The  second  detection  system  is  a  fluorescer  detection  system  and  produces  an 
open  response  function.  This  second  set  of  response  functions  allows  the  experi¬ 
menter  to  determine  the  measured  signals  using  a  different  set  of  response  func¬ 
tions  and  thus  reduce  the  possibility  of  errors  in  the  measured  signals  from  the 


response  functions.  The  inband  response  function  for  this  system  is  defined  to  be 
of  the  same  width  as  the  corresponding  closed  response  function.  This  response 
function  is  depicted  in  Figure  2. 


The  next  step  was  to  define  the  open  and  closed  response  functions.  The 
derivation  of  the  response  functions  can  be  found  in  Appendix  A,  and  is  based  on 
a  simplification  of  the  detector  response  function  presented  by  G.M.  Gorbachenko 
and  others  in  reference  (4).  The  response  functions  used  in  the  study  are  asym¬ 
metric  and  present  a  simple,  but  realistic  and  analytically  representable  shape. 
Additionally,  the  response  functions  allow  for  an  unbiased  unfolding  using  various 
basis  functions  to  define  the  unfolded  spectrum  since  ;  the  response  functions  are 
not  modeled  by  planckians  or  cubic  splines.  Equations  (3)  and  (4)  represent  the 
closed  and  open  response  functions  used  in  this  study. 
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where 


R,c  =  the  closed  function  of  the  ith  detector 
R,°  =  the  open  response  function  of  the  ith  detector 
E,°  -*  the  k— edge  of  the  fluorescer  for  the  ith  detector 
E*  ■*  the  k— edge  of  the  filter  for  the  ith  detector 
E  *  energy 
«  ™  the  ith  detector 


With  each  detector  having  a  resolving  power  of  about  1.5,  the  experimenter 
can  achieve  a  fifty  percent  overlap  in  the  inband  response  functions  of  the  detec¬ 
tors  by  using  thirteen  closed  response  detectors.  The  other  seven  detectors  can 
then  be  used  as  open  response  detectors  covering  the  entire  energy  range.  Wjth 
the  detectors  arranged  in  this  manner,  the  experimenter  is  able  to  detect  the  entire 
spectrum  using  the  inband  response  functions  even  if  a  detector  fails.  This  fact  is 
extremely  important  when  fielding  experiments.  For  this  detector  arrangement, 
the  k-edges  of  the  filter  and/or  fluorescer  for  the  detectors  are  as  listed  in  Table  I. 
Thus,  Equation  (3)  corresponds  to  detectors  1  to  13  and  Equation  (4)  corresponds 
to  detectors  14  to  20. 


a.  E,°  and  E,1  are  computed  with  respect  to  E,0 

Definition  of  Spectra 

During  the  radiation  detection  and  deconvolution  process,  three  spectra  are 
required.  The  first  spectrum  is  the  actual  spectrum  (Sa(E)).  The  actual  spectrum 
is  unknown  and  the  goal  of  the  deconvolution  process  is  to  recover  an  approxima¬ 
tion  to  the  actual  spectrum  from  the  measured  signals  as  discussed  below. 
Secondly,  the  process  requires  a  predicted  spectrum  (5p(E)).  This  predicted  spec¬ 
trum  is  based  on  basic  physics  and  on  all  the  available  knowledge  concerning  the 

.1^ 

spectral  shape.  This  predicted  spectrum  is  used  to  calculate  the  measured-to- 


predicted  and  unfolded-to-predicted  ratios.  Thus,  a  wrong  predicted  spectrum 
will  still  allow  a  proper  deconvolution  of  the  approximation  to  the  actual  spec¬ 
trum;  but  in  the  iterative  scheme  and  the  multiplier  spline  approach,  the  detail  of 
the  predicted  spectrum  is  intentionally  retained  (6).  Also,  an  error  could  be  intro¬ 
duced  by  dividing  a  very  small  predicted  signal,  or  a  predicted  signal  of  0.0,  into  a 
large  measured  signal  or  vice  versa.  So  the  best  possible  predicted  spectrum  is 
desired. 

The  final  spectrum  required  is  the  unfolded  spectrum  (S„  ( E )).  This  spec¬ 
trum  is  the  approximation  to  the  actual  spectrum  and  is  produced  during  the 
deconvolution  process  (i.e.  the  unfolding  of  the  Fredholm  equation  based  on  the 
actual  spectrum  as  implicitly  defined  in  Equation  (1)).  This  unfolded  spectrum  is 


given  by: 


where 


S.(£)-  2° ,F,(E,P ) 


nj  =  the  number  of  basis  functions  be  used 
a;  =  the  coefficient  for  the  jth  basis  function 
F j{E ,P)  =  the  jth  basis  function  and  is  a  function  of 
energy  and  some  other  parameter(s)  P 

The  basis  functions  are  the  building  blocks  used  to  construct  the  unfolded  spec¬ 
trum.  This  study  examined  planckians  and  cubic  splines  as  sets  of  basis  func¬ 
tions.  The  planckian  function  is  given  by  Equation  (6)  (1:5.5)  and  the  cubic 
spline  basis  functions  are  four  point  Lagrangian  interpolating  splines.  The  splines 
are  discussed  in  the  next  section  and  are  derived  in  Appendix  B. 


E 


T  —  the  temperature  of  the  black  body  distribution 
k  =“  the  Boltzman  constant 

Flux  Non-Negativity  Constraint 

In  the  detection  of  radiation  spectra,  a  negative  radiation  intensity  at  any 
energy  level  is  physically  meaningless.  However  by  choosing  an  inappropriate  set 
of  basis  functions  for  the  unfolded  spectrum  or  having  faulty  data,  a  negative 
spectral  value  may  be  produced  during  deconvolution.  In  order  to  ensure  a  posi¬ 
tive  spectrum,  a  non-negativity  option  was  included  in  the  computer  program. 
Equation  (7)  represents  the  technique  used  to  ensure  the  spectra  remained  posi¬ 
tive. 

0.0  S„(£)<0.0 

5b,(£)  =  {5U(£)  S„(£)>0.0  (7 

where  Su'  (E)  is  the  constrained  unfolded  spectrum 

Definition  of  Signals 

Using  the  three  spectra  defined  previously,  four  signals  are  defined.  The  first 
is  the  ideal  signal  ( Y*(E).  However,  as  discussed  earlier,  the  ideal  signal  is  unk¬ 
nown  due  to  calibration  error  in  the  response  functions  and  measurement  and 
detection  uncertainty.  These  ideal  signals  are  approximated  by  the  measured  sig¬ 
nals  from  the  detectors  used  during  the  experiment  (  y',m(E)).  These  measured 
signals  contain  the  errors  and  uncertainty  discussed  above.  This  error  and  uncer¬ 
tainty  is  represented  by  <7,,  the  estimated  standard  deviation  of  the  measured-to- 
predicted  ratios  caused  by  the  error  distribution  in  Y,rn.  The  predicted  signal 
(Y,P(E))  is  then  given  by: 

OO 

r.’(E)-  {S,(E)R,(E)dE 


(6 


Rt{E)  »  the  calibrated  response  function  of  the  ith 

detector,  either  open  or  closed  depending  on  the 
detector 

The  final  signal  is  the  unfolded  signal,  ( Y*(E)),  and  is  given  by: 

OO 

Y,-(E)-  jS,(E)R,(E)JE 


Formulation  of  x2  Teat  Statistic 

For  convenience,  the  measurements  in  this  study  are  specified  as  ratios  to  the 
predicted  signals  rather  than  in  engineering  units  (which  could  not  be  interpreted 
without  detailed  knowledge  of  the  experiment).  Thus,  two  ratios  are  defined. 

The  first  is  b, ,  the  measured-to-predicted  ratio.  This  ratio  approximates  the 
ideal-to-predicted  ratio.  The  measured-to-predicted  ratio  is  given  by: 


The  second  ratio  is  the  unfolded-to-predicted  ratio,  c, ,  and  is  given  by: 


Thus,  the  objective  of  the  unfolding  or  deconvolution  process  is  to  choose  the 
a;  ’s  and  parameter(s),  P,  in  the  F;(E,P)'s  of  Equation  (5)  to  minimize  the 
difference  between  c,  (which  depends  on  them)  and  b,  (which  is  dependent  on  the 
measurement  data).  In  order  to  give  appropriately  increased  weight  to  the  more 
accurate  detectors,  a  test  statistic  was  formulated  and  minimized.  Equation 
(12)  defines  the  x2  test  statistic. 


where 


n i  *  the  number  of  detectors 

cr,  =  the  standard  deviation  in  the  ith  measured— to— predicted  ratio 

However,  before  x2  can  be  minimized,  the  above  mentioned  signal  integrals 
must  be  evaluated.  In  this  study,  the  composite  midpoint  rule  was  used  to 
numerically  evaluate  all  integrals.  Since  the  detectors  only  detected  the  spectrum 
above  E j°,  the  integrals  were  evaluated  beginning  at  this  energy.  Also,  the 
integral  was  truncated  at  128*£j°.  The  approximation  to  the  signal  integrals  is 
given  in  Equation  (13). 

128£? 

/  Stt  (E)R,(E)dE  ^  ZSu(Ek)R,(Ek)±Ek  (13) 

£?  *-l 

where 

nk  =  the  number  of  energy  bins 
Ek  —  the  midpoint  energy  for  the  kth  energy  bin 
i\Ek  —  the  energy  bin  width  of  the  kth  bin 

By  selecting  narrow  energy  bin  widths,  the  error  introduced  in  the  unfolding  pro¬ 
cess  because  of  Equation  (13)  becomes  negligible  when  compared  to  the  other 
errors.  Table  II  presents  a  study  using  one  planckian  basis  function  with  a;  equal 
to  2.0  and  a  temperature  of  S.OEf*  to  evaluate  the  unfolded  signal.  Due  to  the 
equal  resolution  of  the  detectors,  only  one  detector  was  used  as  an  example. 

Minimization  of  \2  Teat  Statistic 

Once  x2  was  defined  using  a  trial  spectrum,  an  iterative  process  to  minimize 
\2  was  required.  Three  methods  were  studied.  The  first  method  conducted  the 
minimization  by  a  least  squares  analysis.  This  method  was  useful  for  linear 
optimization;  but,  the  method  was  not  applicable  for  adjusting  the  temperatures 
in  the  planckians  or  adjusting  the  knots  for  the  cubic  splines.  Also,  the 


method  could  not  be  used  with  the  flux  non-negativity  constraint.  Finally,  a  large 
amount  of  computer  time  would  be  used  calculating  inverse  matrices.  The  calcu¬ 
lation  of  the  inverse  of  a  matrix  also  presented  a  possible  problem.  The  matrix 
could  be  singular  or  near  singular.  However,  this  last  problem  could  be  corrected 
by  selecting  the  appropriate  basis  functions.  Based  on  the  need  for  a  flux  non¬ 
negativity  constraint,  the  method  was  rejected. 

The  second  technique  studied  was  the  steepest  descent  method  of  minimiza¬ 
tion.  This  technique  proved  reliable  but  required  a  large  number  of  iterations 
when  the  number  of  parameters  in  the  basis  function  or  the  number  of  basis  func¬ 
tions  was  increased. 

The  Fletcher-Powell  minimization  technique  was  the  third  method  tested  and 
was  selected  for  use  in  this  study.  The  method  consists  of  calculating  the  gradient 
of  the  function,  x2-  and  then  multiplying  this  gradient  matrix  by  a  correction 
matrix  which  modifies  the  search  direction.  A  detailed  description  of  the  method 
can  be  found  in  references  (3)  and  (5:75-78). 

Two  problems  do  exist  with  this  method.  The  Fletcher-Powell  method  was 
designed  for  those  functions  for  which  the  gradient  can  be  determined  analyti¬ 
cally.  However  in  this  study,  the  gradient  must  be  calculated  numerically.  Thus, 
an  additional  error  is  added  into  the  calculation.  Secondly,  the  method  is  used  to 


find  a  local  minimum  and  does  not  guarantee  this  local  minimum  to  be  the  global 
minimum.  However,  during  this  study  the  local  minimum  problem  was  not 
encountered  since  the  basis  functions  used  contained  only  one  local  minimum  or 
global  minimum. 

Once  the  search  direction  was  determined  using  the  Fletcher-Powell  method, 
the  distance  to  move  in  that  direction  had  to  be  calculated.  This  distance  was 
calculated  by  using  a  line  search  routine  to  minimize  \*(param.+t*s).  In  this 
study,  t  is  the  step  size  or  distance  to  move  in  the  search  direction  and  s  is  the 
unit  search  direction  calculated  in  the  Fletcher-Powell  minir  ization  technique. 
Note  that  a  search  direction  is  based  on  the  gradient  of  \ 2  and  is  a  function  of  all 
the  parameters  used  in  the  basis  functions. 

Two  line  search  routines  were  evaluated  during  this  study.  The  first  method 
consisted  of  calculating  the  value  of  \2  at  two  locations  and  comparing  the  values. 
The  lower  value  for  \ 2  was  retained  and  the  search  interval  was  expanded  or 
reduced  in  order  to  further  minimize  x2.  This  method  was  continued  until  the 
interval  between  two  locations  was  less  than  a  given  value.  Finally,  the  total  dis¬ 
tance  travelled  was  calculated.  This  value  was  then  used  as  the  distance  to  travel 
in  the  search  direction. 

The  second  method  evaluated  and  the  one  selected  for  this  study  used  the 
value  and  slope  of  \2  at  two  locations  to  construct  a  cubic  fit.  This  cubic  fit  was 
then  used  to  determine  the  distance  to  travel  in  the  search  direction.  Reference 
(5:76-80)  presents  the  method  in  detail.  One  modification  to  the  method  was 
required.  If  the  functional  values  of  x2  are  equal  and  the  slope  of  the  second  or 
new  location  in  the  search  direction  is  0,  then  the  search  distance  is  set  equal  to 
the  distance  traveled  between  the  original  \2  location  and  this  final  \2  location. 
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Normal  Distribution  for  Simulating  Measurement  Error 


Once  the  basic  unfolding  technique  was  determined,  the  ability  to  add  simu¬ 
lated  measurement  error  to  6,  was  added.  The  method  employed  to  calculate  this 
normal  distribution  for  the  simulated  measurement  error  is  found  in  reference 
(7:949-953)  and  is  given  by: 


where 


IP  5 

n 

V2 


(N) 


Uj  =*  the  jth  number  in  a  sequence  of  psuedo— random 
numbers  distributed  in  the  interval  [0,1] 
n  =*  the  number  of  psuedo— random  numbers  utilized 

For  simplicity,  twelve  values  of  U }  were  used.  These  values  were  obtained  using  a 
random  number  generator  available  on  the  UNIX  computer  system.  A  sample  of 
ten  normal  variates  generated  in  this  way  was  tested  for  skewness  and  kurtosis. 

In  both  tests,  the  distribution  could  not  be  rejected  as  normal. 
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IV.  Development  of  the  Computer  Program 

Introduction 

In  Section  T H,  the  basic  theory  used  to  develop  the  methodology  involved  in 
the  deconvolution  or  unfolding  technique  was  discussed.  This  section  presents  the 
procedure  used  to  convert  this  methodology  into  a  general  computer  program. 

The  program  was  developed  in  five  basic  steps:  development  of  the  test  spectra; 
development  of  a  planckian  unfolding  program;  addition  of  simulated  measure¬ 
ment  error;  development  of  the  cubic  spline  basis  functions;  and  formalization  of 
the  final  program.  Each  step  is  presented  in  this  section  along  with  a  discussion 
of  the  problems  encountered. 

Development  of  Test  Spectra 

In  the  field  experiment  the  signals  are  produced  and  measured  during  the 
experiment;  thus  the  predicted  spectrum.  Sp(E),  the  calibrated  response  functions. 
R,{E),  the  measured-to-predicted  ratios.  bt,  and  the  estimated  standard  devia¬ 
tions  of  the  measured-to-predicted  ratios.  <r, ,  are  known.  However  in  order  to 
conduct  this  study  these  values  had  to  be  simulated.  By  selecting  parameters  for 
a  set  of  basis  functions,  an  actual  spectrum  was  simulated  using  the  same  method 
as  used  to  produce  an  unfolded  spectrum  in  Equation  (5).  The  actual  spectrum 
was  then  folded  with  the  calibrated  response  functions  derived  in  Section  III  to 
simulate  the  measured  signals.  The  predicted  spectrum  and  signals  were  also  pro¬ 
duced  in  this  manner.  Finally,  the  estimated  standard  deviations  for  the 
measured-to-predicted  ratios  were  defined.  This  program  assumes  the  same  stan¬ 
dard  deviation  for  each  detector. 

Planckian  Unfolding  Code 

The  second  step  of  the  computer  development  consisted  of  writing  a  program 
based  on  the  methodology  in  Section  III  that  could  unfold  the  actual  spectrum 


using  the  same  type  of  basis  function  for  the  unfolded  spectrum  as  was  used  to 
construct  the  actual  spectrum.  In  order  to  simplify  the  calculation  required,  the 
integral  used  to  calculate  the  predicted  signal  in  Equation  (8)  was  normalized. 

The  calibrated  response  functions  were  then  multiplied  by  this  normalization  con¬ 
stant.  Thus,  the  measured-to-predicted  and  unfo!ded-to- predicted  ratios  are  sim¬ 
ply  the  measured  signal  and  unfolded  signal.  However,  the  calibrated  response 
functions  that  were  multiplied  by  the  normalization  constant  must  be  used  to  cal¬ 
culate  these  signals. 

To  accomplish  this,  a  set  of  planckian  basis  functions  were  selected  as  the 
basis  functions  for  the  actual,  predicted,  and  unfolded  spectra.  The  computer 
program  used  the  Fletcher-Powell  minimization  technique  to  vary  the  parameters 
and  coefficients  in  the  basis  functions  of  Equation  (5)  and  determine  the  best 
parameters  for  the  unfolded  spectrum.  The  program  was  validated  using  one,  two. 
and  three  basis  functions.  In  all  cases,  the  predicted  spectrum  was  set  equal  to 
the  actual  spectrum. 

This  validation  is  discussed  in  Section  V.  During  this  validation,  the  two 
methods  for  determining  the  search  distance  is  discussed  in  Section  III  were  com¬ 
pared.  Since  the  cubic  fit  was  faster  and  required  fewer  iterations  of  the  line 
search  subroutine,  it  was  selected. 

During  this  validation  a  possible  problem  was  uncovered  even  though  the 
problem  did  not  affect  the  validation.  The  operator  of  the  computer  program 
should  note  that  it  is  possible  for  the  delta  used  in  the  calculation  ('if  the  gradient 
during  the  minimization  of  y2  to  be  greater  than  the  step  size  or  distance  traveled 
in  the  search  direction.  When  this  occurs,  a  passible  error  in  the  calculation  of  \2 
oould  result. 

The  validation  of  the  Fletcher-Powell  method  also  included  a  comparison  of 
the  number  of  iterations  required  before  convergence  for  this  method  and  the 


steepest  descent  method.  The  Fletcher-Powell  technique  was  dramatically  more 
efficient  and  required  fewer  iterations,  especially  with  the  larger  number  of  param¬ 
eters. 

Based  on  this  first  program  and  its  modification,  the  technique  selected  to 
minimize  x2  was  the  Fletcher-Powell  minimization  technique  with  the  modified 
cubic  line  search  routine.  In  addition  the  flux  non-negativity  constraint  discussed 
in  Section  III  was  used  in  this  program. 

Simulated  Measurement  Error 

The  third  step  of  the  program  development  was  the  user  option  to  add  simu¬ 
lated  measurement  error  to  6,  .  The  simulated  measurement  error  was  used  to 
account  for  the  calibration  errors  and  the  errors  from  the  measurement,  recording, 
and  transmission  of  the  ideal  spectrum.  This  was  accomplished  by  inserting  the 
computer  code  necessary  to  calculate  Equation  (14).  The  simulated  measurement 
error  option  was  included  in  the  program  immediately  after  6,  was  calculated 
and  <7,  was  defined. 

Cubic  Spline  Basis  Function 

Once  the  general  program  was  developed  and  validated  using  the  planckian 
basis  functions,  another  basis  function  subroutine  was  added  to  allow  for  the 
actual,  predicted,  and  unfolded  spectra  to  be  constructed  from  either  planckians 
or  cubic  splines. 

The  cubic  spline  basis  function  consisted  of  two  linear  segments,  a  planckian 
tail,  and  a  variable  number  of  cubic  spline  segments.  The  segments  were  defined 
between  two  consecutive  knot  locations.  The  cubic  splines  used  were  cubic 
Lagrangian  interpolating  splines.  The  basis  functions  were  then  defined  as  the 
combination  of  four  segment  functions.  The  details  of  this  cubic  spline  basis  func 
tion  is  presented  in  Appendix  B. 


Two  forms  of  this  new  subroutine  were  considered  and  developed.  The  first 
form  required  the  knots  used  in  determining  the  cubic  splines  to  be  fixed.  This 
form  was  validated  as  shown  in  Section  V.  However,  in  order  to  effectively 
minimize  x2,  the  knots  must  be  considered  variable.  In  other  words,  the  knots 
were  considered  parameters  of  the  cubic  splines.  The  final  version  of  the  program 
used  the  cubic  splines  with  variable  knots  as  a  possible  type  of  basis  function. 

Decon7.f 

Decon7.f  is  the  final  version  of  the  deconvolution  or  unfolding  program 
presented  in  this  study.  The  computer  program  allows  for  planckian  basis  func¬ 
tions,  cubic  spline  basis  functions,  or  other  spectra  (input  from  a  file)  to  be  used 
as  the  actual  and  predicted  spectra  and  either  the  planckians  or  the  cubic  splines 
can  be  used  to  define  the  unfolded  spectrum.  The  program  then  uses  the 
Fletcher-Powell  minimization  technique  combined  with  the  modified  cubic  line 
search  routine  to  minimize  y2  by  modifying  the  parameters  in  the  basis  functions 
used  to  produce  the  final  unfolded  spectrum.  The  next  section  will  present  the 
validation  and  results  of  this  program.  In  addition,  the  documentation  and 
pseudo-code  for  the  program  are  presented  in  Appendix  C  and  the  source  code  is 
presented  in  Appendix  D. 
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V.  Results  and  Discussion 

Introduction 

Section  IV  discussed  the  development  of  the  computer  program,  during  which 
several  validation  cases  were  considered.  This  section  presents  those  validation 
cases  as  well  as  the  results  from  the  preliminary  study  on  the  effects  of  simulated 
measurement  error  in  the  measured  to  predicted  ratios  and  the  weaknesses  of  the 
two  types  of  basis  functions.  For  this  study  a  final  x2  value  of  less  than  or  equal 
to  the  number  of  detectors  used  was  considered  acceptable.  Also,  in  this  section, 
the  initial  spectrum  refers  to  the  initial  guess  at  the  unfolded  spectrum,  and  the 
actual  and  predicted  spectra  are  identical. 

Validation  Cases 

The  validation  cases  presented  were  developed  to  ensure  the  computer  pro¬ 
gram  was  functioning  properly  and  to  test  the  methodology  used  to  construct  the 
computer  program.  The  goal  of  the  validation  cases  was  to  unfold  the  exact  spec¬ 
trum  that  was  used  to  construct  the  measured  signals.  For  these  cases,  an  arbi¬ 
trary  spectrum  was  selected  as  the  actual  spectrum  and  this  spectrum  was  then 
folded  with  the  calibrated  response  functions  to  define  the  measured  signals.  .As 
mentioned  above,  the  predicted  spectrum  was  identical  to  the  actual  spectrum. 
Also,  simulated  measurement  error  was  not  added  to  the  measured-to-predicted 
ratios.  Finally,  another  arbitrary  spectrum  was  selected  as  the  initial  guess  at  the 
unfolded  spectrum.  In  order  to  validate  the  computer  program  and  the  methodol¬ 
ogy,  the  unfolded  spectrum  should  converge  to  the  actual  spectrum. 

Case  BP.  These  cases  represent  the  benchmark  or  validation  cases  for  the 
planckian  basis  functions.  The  three  cases  presented  validate  the  use  of  one,  two, 
and  three  basis  functions.  The  results  of  these  validation  ca-i  -  ire  presented  in 
Table  III.  In  addition,  a  more  detailed  validation  of  the  one  basis  function  case  is 
presented  in  Appendix  E.  In  all  the  planckian  benchmark  cases,  the  unfolded 
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spectrum  was  indistinguishable  from  the  actual  spectrum  when  plotted  as  noted  in 
Table  III  since  the  parameters  of  the  unfolded  spectrum  converged  to  those  of  the 
actual  spectrum. 

TABLE  in 


Validation  Cases  for  Planckian  Basis  Functions  (a,b) 


Case 

Spectrum 

Parameters(c) 

\  / 

X2 

BF1 

BF2 

BF3 

a 

T 

a 

T 

a 

T 

BPl 

Actual 

2.0 

5.0 

Initial 

1.0 

6.0 

LTnfolded 

2.0 

5.0 

9.8E-4 

BP2 

Actual 

3.0 

6.0 

8.0 

10.0 

Initial 

5.0 

2.0 

1.0 

7.0 

Unfolded 

3.0 

6.0 

8.0 

10.0 

5.9E-3 

BP3 

Actual 

8.0 

5.0 

2.0 

10.0 

1.0 

15.0 

Initial 

3.0 

2.0 

-4.0 

6.0 

2.0 

8.0 

Unfolded 

8.0 

5.0 

1.9 

9.9 

1.1 

15.0 

2.6E-2 

a.  Convergence  criteria:  0.01  or  less  than  a  \%  change  in  \2  for  successive 

iterations 

b.  a,  —  .01 

c.  BF  stands  for  basis  function 

Case  BF.  The  first  four  BF  cases  represent  using  cubic  splines  with  fixed 
knots  as  the  basis  function  and  using  the  same  knots  for  the  actual,  predicted, 
and  unfolded  spectra.  These  cases  validated  the  use  of  the  cubic  spline  basis  func¬ 
tions  with  the  methodology  validated  in  the  BP  cases.  The  results  of  these  cases 
are  presented  in  Table  IV.  Once  again,  the  unfolded  and  actual  spectra  were 
identical  when  plotted  and  the  parameters  of  the  spectra  also  converged. 

Four  additional  cases  were  also  considered.  These  last  four  cases  were  used 
to  validate  the  applicability  of  using  fixed  knots.  In  these  cases,  the  knots  used  to 
form  the  actual  and  predicted  spectra  were  different  than  those  used  to  form  the 
unfolded  spectra.  .As  expected,  the  deconvolution  technique  was  unable  to  correct 
for  the  error  in  the  knot  locations. 


Thus,  all  four  cases  yielded  a  x2  of  greater  than  20.  The  table  of  results  and 
plots  for  these  cases  can  be  found  in  Appendix  E.  Based  on  these  results,  the  knot 
locations  were  included  sis  parameters  in  the  basis  functions. 


TABLE  IV 

Validation  Cases  for  Cubic  Spline 
Basis  Functions  with  Fixed  Knots  (a,b,c) 


Case 
BFl  & 

Spectrum 

Actual 

4.0 

Cubic  Coefficient 
5.0  3.0  1.6 

\7’ 

1.0 

/ 

T 

10.0 

X2 

BF2(d) 

BFl 

Initial 

2.0 

7.0 

1.0 

3.0 

4.0 

10.0 

Unfolded 

4.0 

5.0 

3.0 

1.6 

1.0 

10.0 

4.3E-3 

BF2 

Initial 

2.0 

7.0 

1.0 

3.0 

4.0 

20.0 

Unfolded 

4.0 

5.0 

3.0 

1.6 

1.0 

10.0 

4.3E-4 

BF3  & 

Actual 

4.0 

5.0 

3.0 

1.6 

1.0 

10.0 

BF4(e) 

BF3 

Initial 

2.0 

7.0 

1.0 

3.0 

4.0 

20.0 

Unfolded 

4.0 

5.0 

3.0 

1.6 

1.0 

10.0 

2.4E-3 

BF4 

Initial 

2.0 

7.0 

1.0 

3.0 

4.0 

2.0 

Unfolded 

4.0 

5.0 

3.0 

1.6 

1.0 

10.0 

6.6E-4 

a.  Convergence  criteria:  x*<  0.01  or  less  than  a  0.01%  change  in  \2  for  succes¬ 
sive  iterations. 

b.  <7,  =  0.01 

c.  Fixed  knots  at  0.0,1.0,128.0  are  implicit  in  definition  of  these  splines  and  are 
not  listed. 

d.  Knots  for  case  BFl  and  BF2  were  5.0,10.0.25.0.50.0 

e.  Knots  for  case  BF3  and  BF4  were  10.0,25.0,50.0,80.0 

Case  BC.  The  BC  cases  were  formulated  to  validate  the  computer  pro¬ 
gram  using  the  cubic  spline  basis  functions  with  variable  knots.  As  explained  in 
Appendix  B,  these  basis  functions  were  formed  by  combining  four,  four  point 
Lagrangian  interpolating  polynomials.  Each  basis  function  was  a  function  of  the 
knots  selected  to  construct  the  polynomials  and  the  intensity  at  these  knot  loca¬ 
tions.  The  results  of  these  cases  are  presented  in  Table  V.  The  linear  plots  of  the 
actual  and  unfolded  spectra  for  all  cases  are  located  in  Appendix  E. 

These  validation  cases  presented  the  following  noteworthy  points.  First,  the 
study  presented  the  difficulty  the  computer  program  had  in  varying  the  knot  loca¬ 
tion.  The  program  was  only  able  to  determine  an  acceptable  spectra  in  three  of 
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TABLE  V 

Validation  Cases  for  Cubic  Spline 
Basis  Functions  with  Variable  Knots  (a,b,c) 


Case 

Spectrum 

Parameters 

Knots 

Coefficients 

T 

BCl 

Actual 

5.0 

30.0 

5.0 

3.0 

15.0 

60.0 

2.0 

1.0 

Initial 

2.0 

15.0 

1.0 

4.0 

8.0 

35.0 

6 

0,2.0 

Unfolded 

12.0 

13.0 

4.8 

1.5 

13.0 

28.0 

0.44 

1.5 

BC2 

Actual 

5.0 

10.0 

4.0 

5.0 

10.0 

25.0 

50.0 

3.0 

1.6 

1.0 

Initial 

2.0 

15.0 

3.0 

2.0 

5.0 

35.0 

60.0 

5.0 

1.0 

2.0 

Unfolded 

4.0 

11.0 

4.0 

4.9 

10.0 

28.0 

53.0 

2.7 

1.6 

0.88 

BC3 

Actual 

5.0 

10.0 

4.0 

6.0 

10.0 

25.0 

50.0 

5.0 

3.0 

80.0 

2.0 

1.0 

Initial 

2.0 

15.0 

1.0 

3.0 

4.0 

30.0 

45.0 

5.0 

4.0 

60.0 

3.0 

2.0 

Unfolded 

2.4 

13.0 

3.9 

4.8 

12.0 

37.0 

49.0 

4.2 

2.3 

62.0 

2.0 

1.8 

BC4 

Actual 

5.0 

30.0 

5.0 

3.0 

15.0 

60.0 

2.0 

1.0 

Initial 

4.0 

35.0 

4.0 

4.0 

13.0 

50.0 

3.0 

2.0 

Unfolded 

4.5 

30.0 

5.0 

3.2 

15.0 

60.0 

2.0 

1.0 

0.026 


a.  cr,  =  0.01 

b.  Convergence  criteria:  x2^  0-01  or  less  than  0.1%  change  in  x2  for  successive 
iterations 

c.  Fixed  knots  at  0.0,1.0,128.0  are  implicit  in  definition  of  these  splines  and  are 
not  listed. 

the  four  cases.  The  case  that  was  rejected  used  three  fixed  knots  and  three  vari¬ 
able  knots.  This  case  represented  the  minimum  number  of  knots  possible  in  order 
to  have  at  least  one  section  of  the  spectrum  that  is  constructed  from  four  cubic 
splines.  This  difficulty  is  most  likely  due  to  the  fact  that  the  six  knot  cubic  is 
constrained  to  a  linear  or  planckian  fit  in  all  but  one  section.  Thus,  the  small 
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variations  in  the  knot  locations  cause  only  a  small  change  in  x2  and  the  computer 
program  is  finally  stopped  due  to  such  small  changes  in  x2-  More  computer  time 
could  be  used  to  determine  if  the  value  for  x2  is  finally  reduced  but  a  better  solu¬ 
tion  is  to  select  another  set  of  basis  functions,  add  more  knots,  or  change  the  knot 
locations.  This  is  shown  by  the  accuracy  of  the  unfolded  spectrum  for  case  BC4 
when  the  initial  guess  was  close  to  the  actual  spectrum  as  compared  to  case  BCl. 

Finally,  the  validation  results  for  the  cubic  spline  basis  functions  show  the 
need  for  an  accurate  guess  at  the  knot  locations.  By  optimizing  the  selection  of 
knots  for  the  initial  guess  at  the  unfolded  spectrum,  the  cubic  splines  with  vari¬ 
able  knots  would  resemble  the  cubic  splines  with  fixed  knots  and  a  better  unfolded 
spectrum  should  be  produced.  This  optimization  could  be  conducted  by  unfolding 
signals  produced  using  the  predicted  prior  to  unfolding  the  real  data.  However, 
this  will  not  ensure  the  best  knot  are  selected  since  the  actual  spectrum  may  not 
compare  to  the  predicted.  Thus,  the  dependency  of  the  cubic  spline  deconvolution 
technique  on  the  knot  locations  is  verified. 

Test  Cases 

Once  the  final  computer  program  was  validated,  two  preliminary  studies  were 
conducted.  The  first  was  a  study  to  determine  the  effect  of  simulated  measure¬ 
ment  error  in  the  6,  ’ s  on  the  deconvolution  process.  The  second  study  was  con¬ 
ducted  to  demonstrate  the  degree  to  which  the  unfolded  spectrum  would  approxi¬ 
mate  the  actual  spectrum  if  different  basis  functions  were  used  for  the  respective 
spectra. 

Simulated  Measurement  Error  Study.  As  stated,  this  study  consisted 
of  adding  a  simulated  measurement  error  to  6,  in  order  to  simulate  the  measure¬ 
ment  errors  discussed  is  Sections  II  and  IV  and  which  are  present  in  the  experi¬ 
mental  data.  This  simulated  measurement  error  was  constructed  using  Equation 
(14).  As  noted,  the  normal  distribution  was  scaled  by  <x( ,  the  standard  deviation 
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of  the  measured  to  predicted  ratio.  Ten  cases  with  simulated  measurement  error 
and  one  benchmark  case  without  the  simulated  error  were  tested  for  each  type  of 
basis  function.  A  <r,  of  0.15  was  used  in  all  cases.  This  ax  was  a  relative  error 
and  not  an  absolute  error.  In  other  words,  the  cr,  ’s  were  the  percent  of  error  in 
the  b,  's.  The  benchmark  case  was  labeled  case  0  while  all  other  cases  were  num¬ 
bered  according  to  the  number  used  as  the  seed  for  the  random  number  generator. 
Tables  VI  and  VII  represent  the  results  of  this  study.  Figures  3  thru  6  depict  the 
combined  spectra  for  each  case.  Appe  ndix  F  contains  the  individual  plots  of  the 
actual  and  unfolded  spectra  for  each  case. 

TABLE  VI 
Noise  Study  For 


Case 

Parameters 

\  / 

X2 

TF(c) 

BFl 

BF2 

BF3 

a 

T 

a 

T 

a 

T 

Actual 

8.0 

5.0 

2.0 

10.0 

1.0 

15.0 

11.0 

Initial(d) 

3.0 

2.0 

4.0 

6.0 

2.0 

8.0 

9.0 

NP0 

6.4 

4.8 

2.5 

6.6 

2.1 

13.0 

0.010 

11.0 

NP01 

6.1 

4.6 

3.0 

7.5 

2.1 

11.0 

18.0 

11.2 

NP02 

6.2 

4.8 

2.7 

6.3 

1.8 

15.0 

29.0 

10.7 

NP03 

6.5 

4.9 

2.4 

6.8 

2.1 

13.0 

16.0 

11.0 

NP04 

6.6 

5.0 

2.5 

7.2 

1.8 

13.0 

6.0 

10.9 

NP05 

4.8 

4.2 

3.5 

7.0 

2.1 

12.0 

13.0 

10.4 

NP06 

6.1 

4.7 

2.5 

7.0 

1.9 

13.0 

12.0 

10.5 

NP07 

6.4 

4.8 

3.0 

7.0 

1.8 

13.0 

9.7 

11.2 

NP08 

6.0 

4.9 

2.5 

7.3 

1.8 

13.0 

17.0 

10.3 

NP09 

6.2 

5.0 

3.1 

7.0 

1.4 

13.0 

15.0 

10.7 

NP 10 

6.1 

4.7 

2.6 

6.9 

2.0 

13.0 

18.0 

10.7 

Convergence  criteria:  %?<  0.01  or  less  than  a  0.1%  change  in  x2  for  succes¬ 
sive  iterations 

BF  represents  basis  function 

TF  stands  for  the  total  fluence  from  0.0  to  oo 

Same  initial  spectrum  for  all  cases,  results  represent  unfolded  spectrum 
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Figure  4.  Semi-log  Plot  of  the  {Planckian}  Actual  and 
(Pianckian)  Unfolded  Spectra  Versus  Energy  lor  the  NP  Cases 
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Figure  5.  Linear  Plot  of  the  (Spline)  Actual  and  (Spline) 
Unfolded  Spectra  Versus  Energy  for  the  NC  Cases 


__  TABLE  VII 

Noise  Study  for  Cubic  Spline 
Basis  Function  with  Variable  Knots  (a,b) 

Case  Parameters  T  \2 

Knots  Coefficient 


Actual 

5.0 

10.0 

4.0 

5.0 

3.0 

10.0 

25.0 

50.0 

1.6 

1.0 

Input(c) 

2.0 

15.0 

3.0 

2.0 

5.0 

5.0 

35.0 

60.0 

1.0 

2.0 

NC0 

4.4 

11.0 

4.0 

4.9 

2.7 

9.7 

0.020 

27.0 

60.0 

1.6 

0.68 

NC01 

3.8 

12.0 

3.1 

5.3 

2.2 

8.5 

11.0 

30.0 

60.0 

1.8 

0.66 

NC02 

13.0 

18.0 

5.3 

1.4 

2.0 

11.0 

17.0 

35.0 

61.0 

1.1 

0.57 

NC03 

5.7 

8.7 

3.5 

4.8 

3.8 

7.5 

8.6 

42.0 

68.0 

1.4 

0.55 

NC04 

3.5 

9.7 

3.1 

4.5 

2.8 

9.6 

5.5 

33.0 

60.0 

1.7 

0.62 

NCOS 

3.1 

14.0 

4.2 

5.3 

0.92 

7.3 

9.4 

19.0 

68.0 

1.8 

0.49 

NC  06 

3.6 

10.0 

3.1 

5.0 

2.9 

10.0 

8.8 

33.0 

60.0 

1.4 

0.61 

NC07 

5.7 

10.0 

4.1 

5.0 

3.2 

13.0 

7.2 

33.0 

60.0 

1.9 

0.39 

NCOS 

3.2 

9.7 

2.8 

4.5 

2.5 

10.0 

15.0 

31.0 

60.0 

1.6 

0.58 

NCOS 

5.2 

14.0 

4.0 

4.3 

2.0 

11.0 

13.0 

35.0 

56.0 

1.8 

0.51 

NC  10 

3.4 

9.4 

3.6 

4.8 

3.0 

10.0 

18.0 

32.0 

60.0 

1.5 

0.66 

Convergence 

criteria: 

\2<  0.01 

or 

less  than  O.lCc 

change  in 

y  for  successive 

iterations 

b.  Fixed  knots  at  0.0,1.0,128.0  are  implicit  in  definition  of  these  splines  and  are 
not  listed 

c.  Same  initial  spectrum  for  all  cases,  results  represent  the  unfolded  spectrum 
The  results  of  the  preliminary  study  concerning  simulated  measurement  error 

show  that  measurement  error  in  the  measured  signals  has  little  effect  on  the 
unfolded  spectrum.  In  all  but  one  case,  the  final  \2  was  less  than  the  number  of 
instruments  used,  20.  The  relation  between  the  locations  or  amounts  of  measure¬ 
ment  error  and  was  not  determined  in  this  study.  However,  this  could  be  con¬ 
sidered  in  a  future  study.  In  addition  to  this  result,  the  simulated  measurement 
error  study  on  the  planckian  basis  functions  showed  that  the  temperatures  pro¬ 
duced  in  the  unfolded  spectrum  should  not  be  considered  the  actual  black  body 
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temperatures,  for  multi-temperature  spectra.  However  since  only  one  spectrum 
was  considered  a  more  detailed  study  should  be  conducted.  This  is  shown  by  the 
temperatures  arrived  at  for  the  unfolded  spectrum  during  the  simulated  measure¬ 
ment  error  study.  Also,  since  the  planckians  were  normalized,  the  sum  of  the 
coefficients  approximates  the  total  fluence  for  the  spectru  m.  It  should  be  noted 
that  the  total  fluence  is  accurate  to  about  5°c  in  the  planckian  cases  tested.  The 
measured-to-predicted  ratios  are  accurate  to  15 °c  but  some  of  the  ratios  are  too 
large  while  some  are  too  small  so  the  total  fluence  error  is  less  or  the  accuracy  is 
better  for  the  total  fluence  measurement. 

It  should  also  be  noted  that  the  values  for  \2  were  smaller  for  the  cubic 
spline  basis  functions  than  for  the  planckian  basis  functions.  Because  the  cubic 
splines  are  local  functions,  the  splines  can  fit  the  detail  of  the  spectrum  better 
than  the  global  functions  like  the  planckians.  Thus,  the  \2  should  be  lower. 

Applicability  of  the  Basis  Functions  Study.  During  the  unfolding  of 
real  data,  the  actual  spectrum  is  not  known  so  one  can  not  select  the  set  of  basis 
functions  that  were  used  to  construct  the  actual  spectrum.  This  study  was  used 
to  determine  how  well  a  set  of  basis  functions  could  approximate  an  actual  spec¬ 
trum  constructed  from  a  different  set  of  basis  functions.  Two  actual  spectra  were 
approximated  using  an  ideal  situation  in  which  simulated  measurement  error  was 
not  added  to  the  measured  to  predicted  ratios.  Also  in  both  cases  the  non¬ 
negativity  constraint  was  applied.  The  first  was  an  actual  spectrum  constructed 
from  three  planckian  basis  functions.  This  spectrum  was  then  approximated 
using  cubic  spline  basis  functions  with  six,  seven,  and  eight  knots.  The  results  of 
this  study  are  presented  in  Table  VTII  and  the  best  approximation  to  the  actual 
spectrum  is  presented  in  Figures  7  and  8.  All  plots  of  the  actual  and  unfolded 
spectra  are  presented  in  Appendix  F. 


,  / 
\y 


» 

V 


TABLE  Vin 

Results  of  Fitting  a  Planckian  Spectrum 
with  Cubic  Spline  Basis 
Functions  with  Variable  Knots  (a,b,c) 


Case 

BFl(d) 

Parameters 

BF2 

BF3 

\ 

Actual 

4.0  5.0 

Knots 

3.0  10.0 

Coefficients 

7.0 

15.0 

T 

PCI  Initial 

2.0  15.0 

35.0 

0.30  0.50 

0.20 

0.60 

8.0 

Unfolded 

2.3  3.9 

30.0 

-0.28  -4.8 

0.22 

0.17 

12.0 

71.0 

PC2  Initial 

2.0 

35.0 

15.0 

50.0 

0.30 

0.20 

0.50 

0.10 

0.30 

8.0 

Unfolded 

1.0 

34.0 

12.0 

49.0 

-0.43 

0.12 

-3.5 

-0.23 

0.29 

11.0 

38.0 

PC3  Initial 

2.0  10.0 
15.0  35.0 

50.0 

0.30 

0.30 

0.40 

0.20 

0.50 

0.10 

8.0 

Unfolded 

1.8  2.2 

5.9  26.0 

49.0 

0.097 

0.031 

O  6 

CO 
CO  00 

-0.22 

0.18 

12.0 

16.0 

a.  <r,  —  0.15 

b.  Convergence  criteria:  V<  0.01  or  less  than  O.l^j  change  in  \2  for  successive 
iterations 

c.  Fixed  knots  at  0.0,1.0,128.0  are  implicit  in  definition  of  these  splines  and  are 
not  listed 

d.  BF  stands  for  basis  function 


This  portion  of  the  study  demonstrates  the  strengths  and  weaknesses  of  the 
system  very  well.  The  results  show  how  well  the  computer  program  can  fit  the 
tail  of  the  planckian  spectrum.  This  is  due  to  the  fact  the  fit  between  the  last 
two  knots  in  the  cubic  spline  is  defined  as  a  planckian  distribution.  Secondly,  the 
results  demonstrate  how  the  selection  and  quantity  of  knots  effects  the  unfolded 
spectrum  (i.e.  the  use  of  six  or  seven  knots  is  not  acceptable  but  the  use  of  eight 
knots  is  acceptable.). 

Finally,  case  PC3  depicts  the  artifacts  that  can  be  added  to  the  unfolded 
spectrum  as  a  result  of  the  unfolding.  This  case  is  plotted  in  Figures  7  and  8. 

The  final  point  to  note  is  the  negative  coefficients.  These  negative  coefficients  are 
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acceptable  since  the  non-negativity  constraint  was  imposed.  This  demonstrates 
the  fact  the  computer  treats  the  basis  functions  as  merely  mathematical  functions 
and  tries  to  add  or  subtract  them  as  required. 

The  second  "actual"  spectrum  studied  was  composed  of  cubic  spline  basis 
functions  using  seven  kp-:s.  i'he  spectrum  was  approximated  using  one,  two,  and 
three  planckian  basis  functions  for  the  unfolded  spectrum.  The  results  of  the 
study  are  presented  in  Table  IX  and  the  best  approximation  to  the  actual  spec¬ 
trum  is  depicted  by  Figures  9  and  10.  All  plots  of  actual  and  unfolded  spectra 
verses  energy  are  presented  in  Appendix  F. 

The  results  of  this  study  indicate  that  as  the  number  of  planckian  basis  func¬ 
tions  is  increased,  the  value  for  x2  is  reduced.  From  the  results,  there  appears  to 
be  an  optimum  number  of  knots  that  should  be  used  to  fit  a  given  spectrum. 
However  as  seen  from  cases  CP3  and  CP4,  the  addition  of  another  basis  function 
may  not  improve  the  approximation  of  the  actual  spectrum.  Also  note  the  nega¬ 
tive  temperatures  which  have  no  physical  meaning.  However,  the  computer  code 
is  once  again  treating  the  basis  functions  strictly  as  mathematical  functions  and  in 
this  case  the  overall  function  may  be  negative.  In  order  to  require  a  positive  tem¬ 
perature,  a  simple  restraint  on  the  temperatures  in  the  computer  code  could  be 
inserted,  although  this  would  result  in  a  poorer  agreement  with  the  measurements. 

Figure  of  Merit 

In  order  to  evaluate  how  well  the  actual  spectrum  was  approximated  by  the 
unfolded  spectrum,  a  study  of  eight  figures  of  merit  was  conducted.  This  study  is 
presented  in  Appendix  G.  The  goal  of  the  study  was  to  determine  a  figure  of 
merit  that  correlated  with  x2-  Thus,  an  experimenter  would  have  a  reasonable 
idea  of  the  goodeness  of  fit  for  the  unfolded  spectrum.  However,  the  eight  figures 
of  merit  studied  do  not  appear  to  correlate  well  with  x2-  The  main  reason  for  the 
figures  of  merit  not  correlating  with  x2  was  the  calculation  of  the  functions,  x2 


TABLE  IX 

Results  of  Fitting  A  Cubic  Spline  Spectrum 
With  Planckian  Basis  Functions  (a,b) 

Case  Parameters 

Knots  Coefficients 


Actual(c) 

5.0 

10.0 

25.0 

50.0 

4.0 

5.0 

3.0 

1.8  1.0 

BFl(d) 

BF2 

BF3 

Bf4 

a 

T 

a 

T 

a 

T 

a  T 

CPI 

Initial 

3.0 

2.0 

Unfolded 

67.0 

2.5 

CP2 

Initial 

3.0 

2.0 

4.0 

6.0 

Unfolded 

43.0 

1.6 

76.0 

9.9 

CP3 

Initial 

3.0 

2.0 

4.0 

6.0 

2.0 

8.0 

Unfolded 

39.0 

1.5 

40.0 

7.0 

39.0 

12.0 

CP4 

Initial 

3.0 

2.0 

4.0 

6.0 

2.0 

8.0 

1.0  3.0 

Unfolded 

36.0 

-1.3 

40.0 

7.0 

40.0 

12.0 

38.0  1.4 

a.  cr,  =  0.15 

b.  Convergence  criteria:  x2^  0.01  or  less  than  0.1%  change  in  X 2  for  successive 
iterations 

c.  Fixed  knots  at  0.0,1.0,128.0  are  implicit  in  definition  of  this  spline  and  are 
not  listed  BF  stands  for  basis  function 

was  calculated  based  on  the  difference  of  the  signals  and  thus  was  weighted  more 
heavily  at  the  low  energies  while  the  figures  of  merit  were  based  on  spectral  values 
and  received  equal  weight  since  the  energy  bins  were  equal.  Thus,  a  good  x2 
would  not  correspond  to  a  good  figure  of  merit  because  a  small  difference  in  the 
ratios  for  the  last  detector  would  cause  a  large  error  in  the  figure  of  merits  since 
the  error  would  be  spread  over  energies  from  64 E,°  to  128£,°. 
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VI.  Summary  and  Conclusions 

Summary 

Validation  Cases.  The  purpose  of  the  validation  cases  was  to  ensure  each 
element  of  the  computer  code  was  functioning  properly  and  to  note  any  Saws. 

The  three  validation  cases  presented  do  exactly  that.  The  benchmark  cases  for 
the  planckians  and  the  cubic  spline  basis  function  with  fixed  knots  demonstrate 
the  power  of  this  unfolding  technique.  These  cases  show  how  well  an  actual  spec¬ 
trum  can  be  approximated  if  the  type  of  basis  function  used  to  construct  the 
unfolded  spectrum  is  the  best  possible  basis  function  and  no  measurement  error, 
either  real  or  simulated,  is  present.  In  other  words,  if  an  optimum  basis  function 
can  be  determined,  the  computer  program  is  outstanding. 

The  cubic  spline  validation  cases  also  demonstrate  the  dependence  of  the 
unfolded  spectrum  on  the  initial  guess  at  the  knot  locations.  Even  with  a  poor 
choice  for  the  initial  knot  locations,  the  method  will  attempt  to  unfold  the  actual 
spectrum.  However,  in  general  the  variations  in  the  knots  will  be  so  small  as  to 
require  an  unacceptable  amount  of  computer  time  before  the  actual  spectrum  is 
unfolded.  Thus,  it  is  recommended  that  the  predicted  spectrum  be  used  to  deter¬ 
mine  the  optimum  knots  for  the  initial  guess  at  the  unfolded  spectrum.  However, 
this  does  not  guarantee  the  best  unfold  with  the  real  data  so  the  unfold  may  have 
to  be  repeated  using  different  of  knot  locations 

Test  Cases.  The  test  cases  demonstrated  that  error  in  the  measured  sig¬ 
nals  (simulated  measurement  error  in  the  b,'s  in  this  study)  prevent  the  "exact" 
actual  spectrum  from  being  unfolded.  However,  the  general  shape  of  the  actual 
spectrum  was  unfolded  in  all  cases.  The  test  cases  also  demonstrated  the  artifacts 
introduced  in  the  unfolded  spectrum  when  attempting  to  approximate  an  actual 
spectrum  that  is  not  composed  of  the  same  set  of  basis  functions  as  the  unfolded 
spectrum.  Thus,  an  experimenter  can  not  unfold  the  actual  spectrum  but  only  an 
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approximation  to  the  actual  spectrum  and  this  unfolded  spectrum  will  vary 
depending  on  the  error  in  the  data  and  the  type  of  basis  functions  used.  However, 
the  general  shape  of  the  actual  spectrum  is  returned 

Conclusions 

The  results  obtained  from  the  test  cases  and  validation  cases  demonstrated 
the  following  conclusions. 

1.  The  final  convergence  and  the  rate  of  convergence  of  the  unfold  using 
cubic  splines  is  dependent  on  the  knot  locations. 

2.  The  cubic  spline  basis  functions  or  local  basis  functions  yielded  a  lower 
\2  than  the  planckian  or  global  basis  functions  did  when  used  to  unfold 
data  with  simulated  measurement  error. 

3.  The  total  fluence  of  the  spectrum  is  more  accurate  than  the  spectrum. 

4.  The  shape  of  the  spectrum  is  definable  but  the  exact  spectrum  can  not 
be  unfolded. 

5.  The  temperatures  in  the  multi-temperature  planckians  do  not  represent 
the  actual  temperatures  of  the  planckians. 

6.  The  type  of  basis  functions  used  will  have  a  direct  affect  on  the  details 
one  can  unfold  and  on  the  artifacts  added  to  the  unfold  spectra. 

In  addition,  the  validation  cases  and  test  cases  studied  demonstrate  the  use¬ 
fulness  of  this  data  analysis  method.  The  data  analysis  method  be  used  in  the 
following  wavs. 


1.  During  the  experimental  planning  phase,  the  computer  program  could  be 
used  to  determine  the  best  detector  locations  and  response  functions  in 
order  to  achieve  a  reasonable  unfold;  and  to  optimize  the  number  of 
detectors  or  channels  in  order  to  confirm  the  predicted  spectrum  or  to 
determine  variations  from  the  predicted  spectrum. 
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2.  After  the  ideal  detector  locations  are  selected  and  the  response  functions 
are  determined,  the  computer  program  can  be  used  to  evaluate  the  affect 
of  various  basis  functions  on  the  unfold  and  the  affect  of  measurement 
error  on  the  unfold. 

3.  The  computer  program  can  then  be  used  to  unfold  the  real  data. 

4.  Finally,  an  analysis  of  the  propagation  of  erroi  can  be  conducted  using 
the  computer  program  in  order  to  get  uncertainty  bounds  for  the 
unfolded  spectrum  and  to  determine  the  limitations  of  the  unfold. 

In  order  to  convert  Decon7.f  to  handle  actual  data  as  discussed  above,  four 
read  statements  would  have  to  be  inserted.  These  read  statement  would  read  the 
response  functions,  the  predicted  signals,  the  measured-to-predicted  ratios,  and  the 
standard  deviations  of  the  measured-to  predicted  ratios.  Thus,  the  computer  code 
is  capable  of  simulating  test  procedures  or  with  three  simple  modifications,  pro¬ 
cessing  actual  data. 


VTI.  Recommendations 


Recommendations 

The  following  recommendations  are  made  for  continued  study  in  this  area. 
First,  the  simulated  measurement  error  study  should  be  continued  in  order  to 
determine  if  the  simulated  measurement  error  has  a  more  dramatic  effect  on  any 
particular  detector  or  detectors.  In  addition,  this  study  could  consider  the  possi¬ 
bility  of  one  or  more  detectors  failing  and  see  how  this  failure  effects  the  unfolded 
spectrum. 

Secondly,  the  study  should  be  continued  to  see  how  well  the  planekian  and 
cubic  spline  basis  functions  can  approximate  other  spectra.  This  study  could  also 
be  used  to  determine  if  the  artifacts  produced  correspond  to  the  k-edges  or  energy 
ranges  of  certain  detectors. 

The  third  area  of  study  could  be  a  more  detailed  study  of  the  relation  of  y2 
to  a  given  figure  of  merit.  This  study  could  be  conducted  in  conjunction  with  a 
study  to  determine  the  error  or  uncertainty  in  the  unfolded  spectrum.  This  error 
analysis  should  include  the  methodology  used  in  the  deconvolution  process  as  well 
as  the  error  established  in  measuring  the  actual  spectrum. 

The  fourth  area  that  could  be  studied  is  the  effect  on  the  unfolded  spectrum 
of  varying  the  number  of  detectors  used  and/or  their  locations.  These  effects 
could  be  studied  using  various  basis  functions.  Formulation,  implementation,  and 
evaluation  of  new  and  varied  types  of  basis  functions  is  a  fifth  area  for  future 
work. 

Finally,  a  study  could  be  conducted  to  determine  how  well  the  deconvolution 
process  can  resolve  specific  details  in  the  actual  spectrum.  By  selecting  an  actual 
spectrum  and  then  adding  a  known  function  and  varying  its  amplitude  and 
width,  a  test  procedure  could  be  constructed  to  evaluate  the  resolving  power  of 
the  deconvolution  method. 


Appendix  A:  Derivation  of  the  Response  Functions 


Before  the  methodology  used  to  develop  the  program  could  be  implemented, 
typical  response  functions  for  the  detectors  had  to  be  derived.  In  detecting  pulsed 
x-ray  radiation,  a  common  technique  uses  a  filter-fluorescer  detection  system.  In 
principle,  this  system  uses  the  k-edges  of  the  filter  and  fluorescer  to  determine  the 
inband  response  of  the  detector.  This  type  of  detection  system  was  selected  for 
the  closed  response  function  detectors.  The  open  response  function  detectors  used 
the  same  system  without  the  filter. 

The  equation  used  to  derive  these  response  functions  is  (l5)(4-244): 


Jk(E,Q)dVldE  =  J{E)rk{E)ujk 
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where 


Jk(E  ,0)  =  the  intensity  of  the  K  series  excited  by 

radiation  having  an  energy  ranging  from  E  to 
E+dE 

dQ  =  a  solid  angle 

~k{E)  =  the  coefficient  of  photoabsorption  on  the  k 
shell 

.v' k  =  the  fluorescent  yield 

U'  =  the  probability  ofthe  appearance  of  the 

fluorescence  quantum  having  an  energy E, 
d  =  the  foil  thickness 
H{E)  =  reduction  coefficient  of  the  exciting 

radiation  (  considered  mass  attenuation 
coefficient  ) 
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M ,{E)  =  reduction  coefficient  of  the  ith  fluorescence 
line  (  also  considered  mass  attenuation 
coefficient  ) 

For  the  purpose  of  this  study,  Equation  (15)  was  simplified  as  follows.  First 
it  was  assumed  that  the  experiment  was  arranged  to  consider  only  one  angle  and 
that  this  minimized  /i, (cos(?r— 9)).  Also  note  that  in  the  vicinity  of  the  k-edge, 

f‘[E)  is  much  greater  than  /i, .  Thus,  fJ.,(E )  can  be  approximated  by  Th 

cos  d>) 

second  simplification  was  that  only  one  fluorescence  energy  was  studied.  Thus, 
Equation  (15)  can  be  simplified  to: 


Jk(E)  =  {const  ){rk(E)) 
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Next,  t  is  proportional  to  — ~ —  or  (l:Section  1.8.  page  2)  and  near  the 
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k-edge  of  a  material  r  is  proportional  to  /a.  This  can  be  used  to  show  nt(E)  is 


proportional  to 
l'-rom''s: 
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where  E°  is  the  energy  of  the  k-edge.  Thus  Equation  ( 16) 
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g  =  some  constant 

The  open  response  function  of  the  detector  can  then  he  given  by  Equation 
( 1  >  )  f  1:2  15). 


R°(E)  =  f(d)Jk(E) 

where 


(IS 


t(d)  =  the  efficiency  of  the  detector.  (  assumed 
constant  with  energy  ) 


lijt—’l-'lYJJj  <19 

The  constant  and  t(d)  have  been  droppea  due  to  the  normalization  of  the 
predicted  signal  and  the  subsequent  modification  of  the  response  function.  Also 
note  that  the  energy  must  be  greater  than  the  k-edge  of  the  fluorescer  in  order  for 
the  fluorescence  to  occur. 

The  addition  of  a  filter  to  this  open  response  function  detector  only  serves  to 
attenuate  the  arriving  radiation.  Thus  Equation  (19)  becomes: 


where 

n  =  constant 
m  =  constant 

r  =  ratio  of  mass  attenuation  coefficient  values  at 
k— edge 

E 1  =  k— edge  of  filter 

By  conducting  a  study  of  mass  attenuation  coefficients,  r  was  determined  to 
range  from  -1  to  12.  For  this  study  a  value  of  6  was  selected.  Secondly,  a  resolv¬ 
ing  power  of  about  1.5  could  be  achieved  if  the  El  is  twice  the  value  of  E°  for  the 
detection  system.  This  also  established  the  limits  for  the  inband  response  func¬ 
tions. 

Finally,  a  parameter  study  was  conducted  to  ensure  at  least  TOFo  of  the  sig¬ 
nal  would  be  inband  and  that  this  signal  would  still  be  strong  enough  to  be 
detectable  over  the  other  noise,  such  as  scattered  radiation,  in  the  detection  sys¬ 
tem.  Based  on  the  parameter  study,  the  value  of  the  constants  was  selected  as 
follows:  g  =  3.0;  n  =  2.0.  With  a  value  of  2.0  for  n,  and  Ex  =■  2E°,  m  =  0.25 


Appendix  B:  Derivation  of  the  Cubic  Spline 
Basis  Functions 


The  cubic  splines  used  in  this  study  were  four  point  Lagrangian  interpolating 
splines  and  not  the  common  B-spline.  These  cubic  spline  basis  functions  were  a 
combination  of  four  cubic  splines.  Also,  a  linear  functions  or  a  planckian  function 
were  required  in  certain  segments  of  the  spectrum.  The  lagrangian  interpolating 
splines  were  a  function  of  the  knot  locations  chosen  to  construct  the  splines  and 
the  spectral  intensity  at  these  knots.  In  this  study,  the  intensity  was  referred  to 
as  the  expansion  coefficient  for  the  spline  formed  using  that  knot.  The  study 
required  three  fixed  knots. The  first  knot  is  equal  to  0.0;  the  second  knot  is  equal 
to  Ei  and  the  last  knot  is  equal  to  l28.0*Ei  .  Also,  a  minimum  of  six  knots  is 
required  before  a  full  cubic  section  is  defined  by  the  basis  functions. 

The  number  of  basis  functions  is  equal  to  the  total  number  of  knots  minus 
one.  Each  knot  location  except  the  first  and  last  knot  location  is  used  to  define  a 
specific  basis  function.  Additionally,  the  next  to  last  knot's  spectral  value  is  used 
to  define  two  basis  functions.  This  appendix  will  present  the  methodology  used  to 
develop  the  various  basis  functions. 

The  primary  function  used  to  define  the  basis  function  was  the  cubic  spline 
function  given  by  Equation  (21)('2:S5). 


w  here 


La 


■>  (E-knot,) 
^  }}0{knotk-knot,) 

i 


i-’i: 


E  =  energy 
knotk  =  primary  knot 
knot,  —  other  knots  in  the  function 

In  this  study,  four  functions,  LAk(E)'s.  were  used  to  define  the  basis  function 


for  the  primary  knot  locations.  This  basis  function  covered  a  range  from  the 


energy  of  two  knots  to  the  left  of  the  primary  knot  location  up  to  two  knots  to 
the  right  of  the  primary  knot.  In  other  words  if  the  given  knot  is  knot(x),  the 
energy  range  covered  from  knot(x-2)  to  knot(x+2).  Figure  11  depicts  the  forma¬ 
tion  of  a  cubic  spline  basis  function  using  the  four  LAk(E)  's  described  by  Equa¬ 
tion  (21)  and  expanded  in  Equations  (22)  to  (25). 

// (x-xl)(x-x2)(x-x3) 

1  (x4-.xl)(x4-x2)(x4-x3) 

f(o)  =  (x~ x2)(x— x3)(x— x5) 

1  (x4-x2)(x4-x2)(x4-x5) 

r  (2)  =  •  (J-*31)(x-x5)(x-x6) 

(x4— x3)(x4—  x5)(x4  —  x6) 

,  *  =  J  (x— x5)(x— x6)(x— j 7) 

[(x4— xo)(x4— x6)(x4— x7) 

The  solid  line  portions  of  each  of  the  four  LA  k  functions  depicted  in  Figure 
ll(a  to  d)  are  added  to  form  the  final  cubic  basis  function  depicted  in  Figure  lie. 
A  cubic  spline  basis  function  is  formed  in  the  same  manner  for  each  knot  location 
unless  the  function  is  modified  as  discussed  below. 


x2<x  <x3 
x3<x  <  x4 
x4<x  <  x5 
x5<x  <x6 


Three  modifications  were  made  to  these  cubic  functions.  First,  the  detectors 
are  unable  to  detect  energy  below  Ej3  or  above  12S*£1°  .  Thus,  these  values  were 
not  computed  in  the  basis  functions  used  to  define  the  unfolded  spectrum. 
Secondly,  four  complete  cubics  can  not  be  evaluated  between  the  next  to  last  knot 
and  the  last  knot.  Thus,  a  planekian  fit,  as  given  by  Equation  (6).  was  used  in 
this  section.  Also,  four  complete  cubics  could  not  be  calculated  from  knot  two  to 
knot  three  or  after  the  planekian  modification,  from  two  knots  from  the  end  to 
one  knot  before  the  end.  In  order  to  allow  all  other  sections  to  be  fit  as  cubics 
and  to  correct  these  two  sections,  a  linear  fit  was  used. 


In  order  to  visualize  these  modifications,  Figure  12  depicts  the  five  basis  func¬ 
tions  used  for  the  case  consisting  of  six  knots.  The  basis  functions  cover  the  fol¬ 
lowing  knot  locations:  basis  function  1,  1.0  to  50.0;  basis  function  2,  1.0  to  50.0; 
basis  function  3,  10.0  to  80.0;  basis  function-4,  10.0  to  80.0;  and  basis  function  5, 
80.0  to  128.0.  Note  the  error  in  the  functions  at  knot5.  This  is  due  to  the  step 
size  used  to  produce  the  figure.  By  using  a  smaller  step  size  to  calculate  the  basis 
functions,  a  smooth  set  of  basis  functions  can  be  produced  at  this  last  knot. 

Table  X  presents  the  basis  functions  used  and  the  knots  and  type  of  func¬ 
tions  used  to  define  the  basis  functions. 


TABLE  X 

Knots  and  Functions  Used  to  Define  the 
Cubic  Spline  Basis  Functions(a,b) 


Basis 

Primary 

Knot(c) 

Energy 

Knots 

Type 

Function 

Range(d) 

Used 

Function 

1 

2 

2-3 

2,3 

linear 

3-4 

2, 3, 4, 5 

cubic 

2 

3 

2-3 

2,3 

linear 

3-4 

2, 3, 4.5 

cubic 

4-5 

3, 4, 5, 6 

cubic 

3 

4 

3-4 

2, 3, 4, 5 

cubic 

4-5 

3, 4, 5, 6 

cubic 

5-6 

4, 5, 6, 7 

cubic 

4(e) 

5 

3-4 

2, 3, 4, 5 

cubic 

4-5 

3, 4, 5, 6 

cubic 

5-6 

4, 5, 6, 7 

cubic 

6-7 

5, 6, 7, 8 

cubic 

N-3 

N-2 

(N-2) 

h(N-l) 

1  N-3,  N-2,  N-l,  N 

cubic 

N-3 

h  N-2 

i  N-4, N-3, N-2, N-l 

cubic 

N-4 

h(N-3) 

1  N-5,  N-4,  N-3,  N-2 

cubic 

N-2 

N-l 

(N-l 

hN 

N-l.N 

linear 

(N-2 

1-  N-l] 

1  N-3, N-2, N-l, N 

cubic 

(N-3 

h(N-2) 

(  N-4, N-3, N-2, N-l 

cubic 

N-l 

N 

N-l 

hN 

N-l.N 

linear 

N-2 

h(N-l] 

1  N-3, N-2, N-l, N 

cubic 

N 

N 

N-N+l 

N-l 

planckian 

N  equals  the  number  knots  minus  1 

If  less  than  eight  knots  are  used,  the  basis  functions  must  be  modified  to  en¬ 
sure  only  a  linear  fit  is  used  between  knots  2  and  3  and  knots  N  and  N-l. 

Primary  knot  value  is  equal  to  one  in  all  cases.  All  other  knots  are  equal  to 
zero.  Basis  function  is  then  expanded  using  the  expansion  coefficient  for  the 
primary  knot. 

Energy  Range  is  considered  as  energies  from  knot  to  knot 

Also  used  for  basis  functions  5  to  N-4  with  primary  knot  equal  to  basis  func¬ 
tion  plus  one. 
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Appendix  C:  Computer  Program  Documentation 

Prose 

Decon7.F  is  a  fortran  77  computer  program  that  is  written  off  the  UNIX 
ELXSI  BSD  (4.2)  operating  system.  The  computer  program  is  currently  available 
on  the  ICC  computer  system  at  the  Air  Force  Institute  of  Technology.  Decon7.f 
is  a  computer  program  that  can  be  used  to  approximate  the  actual  radiation  spec¬ 
trum  in  the  Fredholm  integral  equation.  Decon7.f  can  also  test  various  types  of 
basis  functions  for  defining  the  unfolded  spectrum.  With  minor  modification,  the 
response  functions  or  number  of  detectors  utilized  in  the  unfolded  process  can  also 
be  changed.  This  allows  the  program  to  be  used  to  study  the  effect  of  the  detec¬ 
tors  and  the  type  of  response  function  on  the  unfolding  process.  The  methodology 
used  to  develop  this  program  is  discussed  in  detail  in  Section  III. 

The  program  is  based  on  three  main  assumptions.  First,  the  energy  range 
considered  is  from  E  j3  to  128£j°.  Secondly  if  the  actual  and  predicted  are  defined 
by  the  computer  program  vice  being  input,  the  program  assumes  the  calibrated 
and  actual  response  functions  are  equal.  Finally,  decon7.f  assumes  20  detectors 
are  utilized  to  measure  the  spectra. 

Currently,  the  program  is  limited  to  planckian  and  cubic  spline  basis  func¬ 
tions  for  defining  the  unfolded  spectrum.  However,  a  new  basis  function  subrou¬ 
tine  can  easily  be  inserted  for  defining  the  unfolded  spectrum.  .Also,  since  the 
function  is  non-linear,  the  possibility  exists  for  the  computer  program  to  begin  to 
loop  between  certain  values  of  y2-  This  could  be  corrected  by  either  adding  a 
maximum  number  of  iterations  allowed  or  by  stopping  the  computer  program  and 
then  restarting  the  program  with  a  higher  convergence  criteria  for  x2. 

Decon7.f  is  run  by  compiling  and  then  executing  the  source  code.  Decon7  will 
then  prompt  the  user  for  all  input  required  and  list  all  options  available. 
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Pseudo  Code  for  Main  Program 
Input  t  P . 

Input  energy  bin  width  for  calculating  integrals. 

Define  E,°  and  Etl  for  all  detectors  with  respect  to  E P . 

Calculate  response  functions  using  Equations  (3)  and  (-1). 

Input  type  of  predicted  spectrum  to  be  used,  (planckian,  cubic  or  other.) 
If  planckian,  then 

Input  number  of  basis  functions 

Input  coefficient  and  temp,  for  each  basis  function 

Calculate  predicted  spectrum  by  calling  bbfunc 
End  if 

If  cubic  then 

Input  the  number  of  knots 
Input  knot  locations 
Input  expansion  coefficients 
Input  the  planckian  temperature 

Calculate  predicted  spectrum  by  calling  csfunc 
End  if 

If  other  then 

Read  predicted  spectrum  from  input  file 
End  if 


Calculate  Yf  using  Equation  (8) 
Divide  R,  (E )  by  Y,p  to  normalize  Yf 


Define  the  actual  spectrum  (use  the  same  method  as  defining  the  predicted 
method). 

Calculate  6,  using  Equation  (10) 

Input  <7,  for  all  detectors 

Option  to  add  normal  noise  distribution  to  6, .  (Input  1  for  yes) 

If  yes  then 

Input  seed  number  for  random  number  generator 
Calculate  Zx  using  Equation  (14). 

Add  Z,to  6, 

End  if 

Option  for  flux  non-negativity  constraint  (Input  1  for  yes) 

Input  type  of  basis  function  to  define  the  unfolded  spectrum  (planckian  or 
cubic) 

If  planckian  then 

Input  the  number  of  basis  functions 

Input  the  coefficient  and  temperature  for  each  basis  function 
End  if 

If  cubic  then 

Input  the  number  of  knots 
Input  the  knot  locations 
Input  the  expansion  coefficients 
Input  the  planckian  temperature 
End  if 

Initialize  the  H  matrix  for  the  Fletcher-Powell  method  to  the  identity  matrix. 


Input  the  initial  guess  at  the  functional  lower  bound 
Calculate  the  unfolded  spectrum  by  calling  bbfunc  or  csfunc 
Calculate  c,  using  Equation  (11) 

Calculate  \2  using  Equation  (12) 

While  x*  >  convergence  criteria 
Call  minimize 
Recalculate 

Print  iteration  \2,  and  parameters  to  screen 
End  while 

If  \*<  convergence  criteria 

Write  actual  and  unfolded  spectrum  into  an  output  file 
End  if 

End  of  main  program 

Pseudo  Code  for  Subroutine  bbfunc 

Calculate  spectrum  using  equation  (6)  and  apply  flux  non-negativity  option 
if  applicable 

Return  spectrum  values  to  calling  module. 

Pseudo  Code  for  Subroutine  csfunc 

Calculate  spectrum  using  cubic  splines  in  Appendix  B  and  apply  flux  non¬ 
negativity  option  if  applicable. 

Return  spectrum  values  to  calling  module,  flux 


Pseudo  Code  for  Subroutine  minimize  (  5:75-76) 

For  p=l  to  np 

Call  gradient  c, 

Calculate  gradient  of  \2 
Next  p 

Multiply  gradient  matrix  by  H  matrix  to  find  search  direction 
Normalize  search  direction  to  find  unit  search  direction 
Call  subroutine  linesearch 

Reset  II  matrix  to  identity  matrix  if  slope  of  \‘>  0.0  and  recalculate  search 
direction. 

Calculate  new  parameters  for  unfolded  spectra 
Calculate  new  H  matrix 

Call  Lbfunc  or  csfunc  to  calculate  new  spectrum  based  on  new  parameters 
Calculate  pc,  for  new  parameters  using  equation  (11) 

Let  pc,  —c, 

Return  new  parameter  and  c,  to  main  program, 

Pseudo  Code  for  Subroutine  linesearch  (  5:77-80) 

Calculate  functional  value  (ho)  and  slope  (mo)  with  given  parameters 
If  slope  >  0.0 

Return  and  reset  II  matrix 
End  if 

Calculate  t2 

Calculate  functional  value  (h2)  and  slope  (m2)  using  parameter  +  t2  *  unit 
search  direction. 

C-5 


If  m2  >  0.0 


Use  the  two  slopes  and  functional  values  previously  calculated  to  fit  a 
cubic  and  determine  the  line  search  distance 
End  if 

If  m2  <  0.0  and  h2  >  hO 

Cut  the  search  interval  in  half  and  recalculate  m2.  h2 
End  if 

If  m2  =  0.0  and  h2  <  hO 

Line  search  distance  equals  t2. 

End  if 

If  m2  <  0.0  and  h2  <  hO 

Double  search  distance  and  let  mO  =  m2.  hO  =  h2,  tO  =  t2 

Recalculate  m2  and  h2 
End  if 

If  m2  =  0.0  and  hO  =  h2  (modification  to  reference) 

Line  starch  distance  equals  to 

End  if 

Return  line  search  distance  to  minimize  subroutine 

Pseudo  Code  for  Subroutine  calcfunc 

Calculate  new  parameters 

(tparameter  =  parameter  +  t*  unit  search  direction) 

Call  bbfunc  or  csfunc  to  determine  spectrum  values 

Calculate  c,  based  on  new  parameters  using  equation  (11) 

Calculate  \*  using  equation  (12). 


For  p  =  l,np 

Call  gradient  c, 

Calculate  gradient  of  \2 
Next  p 

For  p  =  I,np 

slope  =  slope  +  gradient  x2  *  unit  search  direction 
Next  p 

Return  slope  and  function  value,  to  linesearch  subroutine 

Pseudo  code  for  Subroutine  gradientc 

Let  parameter  =  parameter  +  delta*  parameter  with  respect  to  the  given 
parameter 

Call  Lbfunc  or  csfunc  to  determine  new  spectrum  values 
Calculate  new  c  \  let  it  equal  /,  using  equation  (li) 

Calculate  gradient  c,  using  c,  and  /, 

Return  the  value  for  the  gradient  of  c 

Source  Code  for  Decon7.f 

See  Appendix  D 

List  of  Variables  used  in  the  Source  Code 

See  tab  1  to  this  appendix 

Inputs 

Inputs  required  for  Decon  7.f  are  listed  in  Table  XL 


TABLE  XI 

List  of  Input  Variables 
Var.  Definition 

Ei  k-edge  of  fluorescor  for  1st  detector 

y  energy  bin  width  for  calculating  integrals 

type*  type  of  basis  function  used  to  define  predicted 

spectrum 

duml*  #  of  planckian  basis  functions;  number  of  knots 

for  cubic;  or  stop  feature  for  reading  input 
file 

pparam(j)*  coefficient  and  temperatures  for  planckian 
basis  function;  or  knots,  expansion 
coefficients  and  temperature  for  cubic 
duml  toggle  for  simulated  measurement  error 

dum2  seed  for  random  number  generator 

non  toggle  for  non-negativity  option 

type  type  of  basis  function  used  to  define  the 

unfolded  spectrum 

crt  standard  deviation  for  bt 

nj  number  of  planckian  basis  functions 

nkn  number  of  knots  for  cubic  basis  function 

param(j)  coefficients  and  temperatures  for  each 
planckian  basis  function  or  knots, 
expansion  coefficient  and  temperature  for 
cubic  basis  function 

1  initial  guess  at  the  functional  lower  bound. 

*  also  required  for  the  actual  spectrum 

Outputs 

The  output  occurs  in  two  forms.  First,  the  values  for  \*  and  all  parameters 
are  printed  to  the  computer  screen  during  each  iteration  as  are  the  final  parame¬ 
ters.  Secondly,  the  actual  spectrum  and  unfolded  spectrum  values  at  each  energy 
location  are  printed  into  an  output  file  labeled  output. 


File  Structure 

The  program  uses  one  output  file  and  can  use  two  input  files.  In  order  to 
run  the  program  both  input  files  must  be  defined.  However,  the  two  input  files 
may  be  empty  if  they  are  not  required  to  define  the  actual  or  predicted  spectrum. 
The  file  structure  for  the  output  file  is  a  two  column  table  with  energy  and  inten¬ 
sity  at  that  specific  energy.  The  first  half  of  the  file  contains  the  actual  spectrum. 

These  energies  begin  with  the  maximum  energy  —  and  decrease  in  steps 


of  the  bin  width  (y)  to  a  value  of  E t°  +  -^-.  The  second  half  of  the  file  consists  of 

the  unfolded  spectrum  and  begins  with  energy  E®  +-¥-  and  increases  in  steps  of  y 
to  the  maximum  energy.  The  file  structure  for  both  input  files  consist  of  a  single 
column  of  energy  intensities  for  energies  from  £[°  +-^-  to  the  maximum  energy  in 
steps  of  y,  the  bin  width. 

Validation  Cases 

See  Sections  V  and  VI. 


Test  Cases 


-ee 


V'Ctions  \ 


am 


VI. 


Tab  1:  List  of  Variables  (a) 

Variable  Svmbol  Definition 

Modules!  b) 

in  code 
a 

in  Text 

dummy  matrix  used  to  calculate  H 

4 

aknotl 

-- 

knot  0.0  for  actual  spectrum 

1 

aknot 1 

— 

knot  Eo(l)  for  actual  spectrum 

1 

aknot 15 

— 

knot  128Eo(l)  for  actual  spectrum 

1 

aparamlj ) 

— 

jth  parameter  for  actual  spectrum 

1 

b(i) 

b, 

measured-to-predicted  ratio 

1  .-4.5.6 

c(i) 

c> 

unfolded-to-predicted  ratio 

1.1 

c(i) 

c(i)  for  appropriate  spectrum 

4 

chi2 

o 

V 

chi  squared 

1 

coef(p) 

expansion  coefficient  knot  p-1 

3 

const 

— 

expansion  coefficient  for  planckian 

3 

const 

— 

dummy  variable  for  summing 

1 

d(i) 

dx 

random  numbers 
integral  of  Ri(E)  Sp(E)  dE 

1 

delta 

step  size  to  calculate  gradient 

i 

dum 

— 

dummy  variable 

3 

d  u  m  1 

-- 

dummy  variable 

1 

dum  1 

— 

dummy  variable 

5 

durn  1  ,dum2 

— 

dummy  variables 

1 

d  u  m  3 

— 

%  change  inx2 

1 

dummy 

— 

dummy  variable 

1 

e 

E 

energy  counter 

1.2.3 

eo(i) 

Et° 

k-edge  of  fluorescer  for 

1.2.3. 

e  1  ( i ) 

Etl 

ith  detector 

k-edge  of  filter  for  ith  detector 

4  .5.6,7 

1 

f(p) 

— 

c(i)  for  newpar(p) 

t 

e(p) 

-- 

gradient  of  \*  for  old  parameters 

1 

gprime(  p ! 

-- 

gradient  \*  for  new  parameters 

1 

grade 

-- 

gradient  of  c(i)  for  appropria'e  spectrum 

4 

grade!  p) 

— 

gradient  c(i)  with  respect  to  parameter  p 

1 

grade  p ) 

-- 

gradient  of  tc(p) 

G 

grchi'_’;  p  1 

- 

gradient  of  \*  for  tparam 

6 

grchi'.’lp) 

- 

gradient  \*  with  respect  to  parameter  p 

1 

h 

-- 

correction  matrix  H 

1. 1.3 

ho 

-- 

functional  value  at  (Param-»-to*s) 

3 

ht 

-- 

dummy  functional  value  variable 

3.6 

h2 

— 

functional  value  at  (Pararn  +  t2*s) 

3 

i 

-- 

counter 

1 .1.0.7 

it 

— 

iteration  counter 

1 

j 

-- 

counter 

1 .2,3,4 

k 

— 

counter 

1.2.3. 

kc 

planckian  conversion 

•1.6.7 

1.2.3. 

knf  p) 

constant  (temp,  to  energy) 
knot  number  p 

■1 .3.6.7 

3 

k  rdla) 

— 

dummy  variable  for  knot  location 

4 

knot  1 

knot  0.0.  unfolded  spectrum 

1.3.4. 

knot  2 

knot  Eo(l),  unfolded 

5,6.7 

1.3.4. 

* 

spectrum 

C-l-1 

3.6,7 

Variable 
in  code 

Symbol 
in  Text 

Definition 

Modules(b) 

knotl5 

— 

knot  128Eo(l),  unfolded 
spectrum 

1,3,4, 

5,6,7 

1 

— 

functional  lower  bound 

1,4,5 

m 

- 

dummy  slope  variable 

5,6 

mo 

- 

slope  at  j 

param  +  to*s) 

5 

mu 

— 

variable  ; 

tar  cubic  fit 

5 

m2 

— 

slope  at  ( 

param  -f  t2*s) 

5 

newpar(p) 

- 

paramfp' 

i  +  (delta)*(param(p)) 

7 

ni 

number  of  detectors 

1,4,5, 

6,7 

nj 

— 

number  of  planckian  basis  functions 

1,2,4, 

nk 

— 

number  of  energies 

1,2,3, 

4.5, 6, 7 

nkn 

— 

number  of  knots  unfolded 
spectrum 

1,3,4, 

5,6,7 

np 

— 

number  of  parameters 

1,4,5, 

6,7 

1,2,3 

4.5. 6, 7 

non 

- 

toggle  for  flux  non-negativity 
option 

norms 

— 

norm  of  sprime 

4 

num 

- 

dummy  variable 

3 

o 

- 

counter 

1,4 

P 

— 

counter 

1,4. 6. 7 

parain(j) 

-- 

jth  parameter  for  unfolded 
spectrum 

1,2.3, 

4, 5, 6, 7 

pc(i) 

- 

c(i)  for  new  parameters 

4 

Pi 

7T 

P‘ 

1.2,3. 

4.5. 6. 7 

pknot  1 

- 

knot  0.0  for  predicted  spectrum 

1 

pknot2 

- 

knot  Eo(l)  for  predicted  spectrum 

1 

pknot 15 

- 

knot  128Eo(l)  for  predicted  spectrum 

1 

planck 

- 

planckian  distribution 

3 

pparam(j) 

- 

jth  parameters  for  predicated  spectrum 

1 

pprime(p) 

— 

new  parameter 

4 

q 

— 

dummy  matrix  used  to  calculate  H 

4 

r(i.k) 

R,(E) 

response  function  of  ith 
detector  kth  energy 

1.4.5, 

6.7 

rand 

- 

random  number  generator 

1 

reset 

— 

reset  counter  for  H  matrix 

4,5 

s(p) 

— 

unit  normal  search 
direction  for  parameter  p 

4.5,6 

sa(k) 

Sa(E) 

actual  spectrum  kth  energy 

1 

sig 

- 

dummy  input  for  <7, 

1 

sigma(i) 

standard  deviation  in  6, 

1.4, 5, 6 

slope 

slope  for  linear  fit 

3 

sp(k) 

S~(E ) 

predicted  spectrum  kth  energy 

1 

spriine(p) 

— 

search  direction  for  parameter  p 

4 

srand 

— 

initializes  random  number  generator 

1 

sufk) 

SAE) 

unfolded  spectrum  kth  energy 

1.2,3 

Mil  k ) 

— 

dummy  variable  unfolded  spectrum 

4 

^ 1  i  1  k ) 

- 

unfolded  spectrum  for  tparamfpl 

6 

-  .  kl 

unfolded  spectrum  for  newpar(p) 

C-l-2 

i 

I 


■■■A 


.■V 

» 

v  -j 

V.-J 

•v'; 

► 


a 


*  ?  *  s/ 


Modules(b) 


Variable 

Symbol 

Definition 

Modules( 

in  code 
sumi 

in  Text 

dummy  variable 

4 

sumj 

— 

dummy  variable 

4 

sumk 

— 

dummy  variable 

4 

sumo 

— 

dummy  variable 

4 

sump 

- 

dummy  variable 

4 

sumran 

— 

sum  of  random  numbers 

1 

t 

- 

temperature  of  planckian 

3 

t 

— 

step  size  in  search  direction 

4,5,6 

tc(p) 

— 

c(i)  for  tparam(p) 

6 

to 

— 

step  size  used  to  determine  ho, mo 

5 

tparam(p) 

— 

parameter(p)  +t*s 

6 

type 

— 

type  of  basis  function 

1,4,5, 

6,7 

5 

tl 

_ 

calculated  step  size 

t2 

— 

step  size  used  to  determine  h2,m2 

5 

V 

— 

dummy  matrix  used  to  calculate  H 

4 

w 

— 

variable  for  cubic  fit 

5 

X 

— 

dummy  variable 

1.2,3 

X 

parameter  transfer  counter 

7 

V 

— 

energy  bin  width 

1.2.3. 

z 

dummy  matrix  used  to  calculate  H 

4.5, 6, 7 
4 

z 

— 

variable  for  cubic  fit 

5 

If  a  parameter  is  listed  once  for  several  modules,  the  parameter  is  passed 
between  the  modules.  If  the  parameter  is  listed  separately  for  each  module, 
the  parameter  is  only  common  to  that  module. 

Module  listing:  1,  Main;  2,  bbfunc;  3,  csfunc;  4,  minimize;  5,  linsearch;  6, 
calcfunc;  7,  gradientc 


Appendix  D:  Computer  Program  Source  Code 


THIS  PROGRAM  IS  THB  FINAL  VERSION  OF  DBCON .  IT  USBS 
THB  FLBTCHER-POWBLL  MBTHOD  FOR  MINIMIZING  CHI2  WHILE 
USING  A  CUBIC  LINB  SEARCH  ROUTINE.  ALL  INPUTS  ARB 
RBQUIRBD  TO  BB  INPUT  FROM  THB  KEYBOARD.  THB  CUBIC 
BASIS  FUNCTIONS  INCLUDB  THB  KNOTS  AS  PARAMETERS. 

THIS  PROGRAM  USBS  NUMERICAL  DBCONVOLUTION  TO 
APPROXIMATE  AN  ACTUAL  RADIATION  SPBCTRUM  WHEN  GIVEN  THB 
MBASURBD  SIGNAL,  A  PREDICTED  SPBCTRUM,  AND  RBSPONSB 
FUNCTION  OP  THB  DBTBCTOR(S).  THB  PROGRAM  ASSUMBS 
13  CLOSBD  RBSPONSB  DBTBCTORS  AND  7  OPBN  RBSPONSB 
DBTBCTORS.  THB  PROGRAM  USBS  AN  BNBRGY  RANGB  OP  Bod) 

TO  128*Bo< 1  ) 

real*8  pparam< 30  ),pknotl,pknot2,  pknotl5 , aparamC 30 > 

real*8  aknotl,aknot2,aknotl5 

real *8  chi2,eo( 20  ), el( 20  >,r( 20,1280  >,e,pi ,kc 

real*8  sa<  1280  ),sp<  1280  ),d<  20  >,b(  20  >,sigma(  20  > 

real*8  sumran,slg,dum3, dummy, x,su( 1280  > 

real* 8  l,param( 30  ),c<  20  ),y,h<  30, 30  > 

real *8  knot 1 , knot 2, knot 15 , srand, rand, const 

integer  nkn , type , i t , non , i , k , 3 , n i , nk , n  3 , duml , dum2 , o , p , np 

common  eo<  20  ) 

open  < unit-7, file-'output' , status- 'new*  ) 
open  < unit-2,  file- ' inputl ', status- 'old*  ) 
open  ( unit-3,  file-' input2’ , status- ' old'  > 
rewind  ( 2  ) 
rewind  <  3  ) 

print  *, 'INPUT  BO< 1  )  ' 
read  *,eo( 1  ) 

NI-THB  NUMBER  OP  INSTRUMENTS  OR  DBTBCTORS 
NK-THB  NUMBER  OP  ENERGY  BINS  USED  FOR  INTEGRATION 
Y-THB  WIDTH  OF  THE  ENERGY  BIN 
ni-20 

print  *,  'INPUT  THB  DESIRED  BNBRGY  BIN  WIDTH,  MINIMUM' 
print  *,  'BIN  WIDTH  IS  0.1* 
read  *,y 

nk=int<  (  128.0*eo<  1  )-eo(  1  )  >/y  ) 

pi-3.1415 

kc-1.0 

DBF INB  ALL  BO'S  WITH  RBSPBCT  TO  BO( 1 > 
do  5  i-2,7 

eo< i >-2.0*eo( i-1 ) 
continue 
eo<  8  )-l .  5*eo<  1 ) 
do  10  i-9,13 

eo< i  )-2*eo( i-1  ) 
continue 
do  15  i-14,20 
eo(  1  )-eo<  i-13  ) 
continue 


c 


DEPINB  Bl< I  )WITH  RBSPBCT  TO  EO< I  ) 
do  20  i-l,ni 

el(  i  )-2.0*eo(  i  ) 

20  continue 

c  DBF INB  ALL  RBSPONSB  FUNCTIONS  AT  SPBCIFIBD  BNBRGIBS 

do  30  i-l,ni 
e-eo(  1  >-y/2.0 
do  25  k-l,nk 
e*e+y 

If  (e.gt.eo(l))  then 
if  < i . It . 13 >  then 

if  (e.lt.el(i))  then 

r(  i,k>-(  1.0/e  >*(  1 .0-exp( -2.0*( eo(  i  >/e>**3  ) ) 
r<  i ,k  )-r (  i ,  k  >*exp(  -0.25*(el(  i  )/e  )**3  > 
else 

r(  i,k  >«(  1.0/e  >*<  1.0-exp(  -2.0*(  eo<  i  )/e  )**3  )  ) 
r<  i  ,  k  >*r  (  i,k)*exp(-1.5*(el(  i  )/e  )**3  > 
end  if 
else 

r(  i ,k  )“<  1.0/e  >*(  1.0-exp(  -3.0*(  eo<  i  >/e  )**3  > ) 
end  if 
else 

r<  i,k>“0.0 
end  if 
25  continue 

30  continue 

c  DEFINE  THE  PREDICTED  SPBCTRUM 

print  INPUT  THE  TYPB  OF  FUNCTION  TO  BE  USBD  TO  MODEL 

print  * , *  THE  PREDICTED  SPECTRUM* 

print  1**PLANCKIAN  BLACK  BODIBS* 

print  *,*  2-CUBIC  SPLINBS* 

print  3-OTHBR* 

read  *,type 

if  (type.eq.l)  then 

print  'INPUT  THB  NUMBBR  OF  PLANCKIAN  BASIS* 
print  *, ’FUNCTIONS  TO  USB* 

print  * , * <  MUST  BB  AN  INTBGBR  LBSS  THAN  15)* 
read  *#duml 

print  ‘INPUT  THB  COBF.  AND  TBMP .  FOR  BACH  BASIS* 
print  * , ’FUNCTION  SBPBRATBD  BY  A  COMMA  (TBMP  IS  IN* 
print  *,  ’UNITS  OF  B0<1>>* 
do  33  j-l,duml 

print  *, ’BASIS  FUNCTION  NUMBBR* , j 
read  *,pparara< j  ), pparam< J+durol > 

33  continue 

non-1.0 

call  bbf unc(  ppar am, sp, y , nk, duml ,pi , kc,non > 
end  if 

if  (type.eq.2)  then 

print  *,* INPUT  THB  NUMBBR  OF  KNOTS  TO  BB  USBD  TO  FIT* 
print  *, *THB  CUBIC  SPLINBS  (MUST  BB  AN  INTBGBR* 
print  *, *6<-iKNOTS<-15 >  NOTBsO.O* 

D.2 


print  *,‘IS  COMPUTBD  AS  THB  FIRST  KNOT,  BO< 1  )  IS* 
print  *, ‘COMPUTBD  AS  KNOT  2,  AND  128.0*BO(1>  IS* 
print  *, 'COMPUTBD  AS  THB  LAST  KNOT.  INCLUDB  THBSB * 
print  *, ‘KNOTS  WHBN  COUNTING  YOUR  KNOTS, AND  BNSURB* 
print  * , * THB  KNOTS  ARB  INPUT  FROM  THB  MINIMUM  TO  THB 
print  *,‘  MAXIMUM. • 
read  *, dural 

print  *, ‘INPUT  THB  KNOTS* 
pknotl-0 .0 

print  *, 'KNOT  l-‘,pknotl 
pknot2-eo< 1 ) 

print  *,'KNOT  2-',pknot2 
do  37  J-l ,duml-3 

print  *, ‘input  knot  ’,J+2 
read  *,pparam< J  ) 
continue 

pknotl5-128.0*eo< 1 ) 
print  *, ‘KNOT* ,duml, *-* ,pknotl5 
do  38  3“l,duml-2 
if  ( j .ne. 1  )  then 

print  *, 'INPUT  THB  EXPANSION  COBF .  FOR* 
print  *,'  KNOT- ' ,pparam(  3-1 ) 
else 

print  *,' INPUT  THB  BXPANSION  COBF.  FOR’ 
print  *,‘  KNOT-' ,pknot2 
end  if 

read  * ,pparam< duml-3+ j  ) 
continue 

print  *, ‘INPUT  THB  TBMP.  TO  BB  USBD  TO  CONSTRUCT  A‘ 
print  * , ‘ PLANCKIAN  FIT  BBTWBBN  THB  LAST  2  KNOTS, • 
print  TBMP  IS  IN  UNITS  OF  Bo‘ 

read  * ,pparam< 2*duml-4 > 
non-1 . 0 

call  csfunc<  pparam,sp,duml,nk,y,non,pknotl, 
fc  pknot2,pknotl5,eo,pi ,kc > 
end  if 

if  (type.eq.3>  then 

print  * , ‘ Sp( B  )  WILL  BB  RBAD  FROM  FILB  inputl,  THB  * 
print  *, 'FILE  SHOULD  CONSIST  OF  SPBCTRUM  VALUBS  FROM 
print  *, '< BO< 1 >+BIN  WIDTH/2.0 )• 

print  * , ‘ TO  128.0*BO(1)  IN  STBPS  OF  THB  BIN  WIDTH* 
print  *, 'INPUT  A  1  TO  CONTINUE  OR  A  2  TO  STOP* 
read  *,  dural 
if  (duml.eq.2)  then 
qo  to  151 
end  if 

e-eo< 1  >-y/2.0 
do  44  k-l,nk 
e-e+y 

read  <  2, *  )  sp(  k  > 
continue 
end  i  f 


CALCULATE  THB  NORMALIZATION  CONSTANT  FOR  D< I >  AND 
RBNORMALIZB  THB  RBSPONSB  FUNCTIONS 
do  55  1-1, nl 
do  45  k-l,nk 

d<  1  )-d(  i  )+sp(  k  )*r<  i  ,k  )*y 
continue 
do  50  k-l,nk 

r(  i ,  k  >-r  (  1 ,  k  )/d(  i  ) 
continue 
continue 

DBFINB  THB  ACTUAL  SPBCTRUM 

print  *, 'INPUT  THB  TYPB  OF  FUNCTION  TO  BE  USBD  TO' 

print  *,'MODBL  THB  ACTUAL  SPECTRUM* 

print  1-PLANCKIAN  BLACK  BODIBS* 

print  2 -CUBIC  SPLINES' 

print  3-OTHBR' 

read  *,type 

if  (type.eq.l)  then 

print  *, 'INPUT  THB  NUMBER  OF  PLANCKIAN  BASIS' 
print  *, 'FUNCTIONS  TO  USB' 

print  * , ' <  MUST  BB  AN  INTBGBR  LBSS  THAN  15)' 
read  *,duml 

print  *, 'INPUT  THB  COBF.  AND  TBMP.  FOR  BACH  BASIS' 
print  *, 'FUNCTION  SBPBRATED  BY  A  COMMA  (TBMP  IS  IN' 
print  *, 'UNITS  OF  B0<1>>* 
do  59  J-l,duml 

print  *, 'BASIS  FUNCTION  NUMBBR ' , } 
read  *,aparam( j >,apararo( J+duml > 
continue 
non-1 . 0 

call  bbf unc<  apar am,sa , y ,nk, duml ,pi ,kc,non  ) 
end  if 

if  (type.eq.2)  then 

print  INPUT  THB  NUMBBR  OP  KNOTS  TO  BB  USBD  TO  FIT 

print  * , ' THB  CUBIC  SPLINBS  (MUST  BB  AN  INTBGBR' 

print  *, ’6<-#KNOTS<-15 >  NOTB:0.0' 

print  *#’IS  COMPUTBD  AS  THB  FIRST  KNOT,  BO( 1 >  IS' 

print  *, 'COMPUTED  AS  KNOT  2,  AND  128.0*BO(1>  IS* 

print  *, 'COMPUTBD  AS  THB  LAST  KNOT.' 

print  * , ' INCLUDB  THESE  KNOTS  IN  THB  COUNT  OF  YOUR* 

print  *, 'KNOTS  AND  BNSURB  THB  KNOTS  ARB  INPUT  FROM* 

print  * , ' THB  MINIMUM  TO  THB  MAXIMUM* 

read  *rduml 

print  *, 'INPUT  THB  KNOTS' 
print  *, 'KNOT  1-0.0* 
aknot 1-0.0 
aknot2-eo( 1 ) 

print  *f'KNOT  2-’,aknot2 
do  62  J-l,duml-3 

print  'input  knot  ',)+2 
read  * ,apararo( j  ) 
continue 


aknotl5-128.0*eo< 1 > 
print  *, 'KNOT* ,dural+3, ,aknotl5 
do  63  J«l,duml-2 
if  <  J .ne.l )  then 

print  * , *  INPUT  THE  EXPANSION  COBF.  FOR* 
print  *, ’KNOT-' ,aparam( j-1 ) 
else 

print  *, *  INPUT  THB  EXPANSION  COBF.  FOR* 
print  * , ’ KNOT- * , aknot2 
end  i  f 

read  * , aparam( duml-3+  j  ) 
continue 

print  *, ’INPUT  THB  TBMP.  TO  BE  USED  TO  CONSTRUCT  A* 
print  * , ' PLANCKIAN  FIT  BBTWBBN  THB  LAST  2  KNOTS,* 
print  * , *  TBMP  IS  IN  UNITS  OF  Bo* 
read  *,aparam< duml*2-4 > 
non-1 . 0 

call  csf unc( apar am,sa , duml ,nk, y ,non , 
aknotl ,aknot2,aknotl5,eo,pi ,kc  ) 
end  i  f 

if  (type.eq.3)  then 

print  * , ' Sa<  B )  WILL  BB  RBAD  PROM  PILB  input2,  THB* 
print  * , 'FILB  SHOULD  CONSIST  OF  SPBCTRUM  VALUBS  PROM 
print  * , ’ <  BO< 1  5+BIN  WIDTH/2 . 0  )  * 

print  * , ' TO  128.0*BO<1>  IN  STBPS  OP  THB  BIN  WIDTH’ 
print  *, ’INPUT  A  1  TO  CONTINUB  OR  A  2  TO  STOP* 
read  *,duml 
if  <duml.eq.2)  then 
go  to  151 
end  if 

e-eo<  1  )-y/2 . 0 
do  69  k-l,nk 
e-e+y 

read  (  2, *  )  sa<  k  ) 
continue 
end  if 

DBF INB  A< I  );  RBCALCULATB  D< I  )  AND  CALCULATE  B( I  ) 
do  75  i-l,ni 
a(  i  )-0.0 
d<  i  )-0.0 
b(  i  )-0.0 
do  70  k-l,nk 

a<  i  )«a<  i  )+sa(  k  >*r<  i, k  )*y 
d<  i  >“d<  i  )+sp(  k  >*r<  1  ,k  )*y 
continue 
b<  i  >-a<  i  )/d<  i  ) 
continue 
DBFINB  SIGMA<  I  ) 

print  *, ’INPUT  SIGMA( I  )  ASSUMBS  ALL  SIGMAS  ARB  BQUAL ' 
read  *,sig 
do  80  i-l,ni 
sigma( i >-sig 


continue 

print  *,'DO  YOU  WISH  TO  APPLY  A  NORMAL/GAUSSAIN  NOISE* 
print  *, ‘DISTRIBUTION  TO  B(I>,  INPUT  A  1  FOR  YES* 
read  *,duml 
if  (duml.eq.l)  then 

print  * , *  INPUT  THE  SBED  FOR  THE  RANDOM  NUMBER' 
print  *, 'GBNBRATOR,  MUST  BE  AN  INTEGBR  * 
read  *,dum2 
const-srand(  dum2  > 
do  82  i-l,ni 
sumran-0 . 0 
do  81  j-1,12 
const-r and(  x  ) 
sumr an -sumr an+const 
continue 

z(  i  )-sigma(  i  >*(  sumr  an- 6 . 0  ) 
b<  i  >-b(  i  >+z(  i  ) 
print  *, 'noise- ' ,z( i > 
print  *, '1*' ,i, 'suraran*' ,sumran 
continue 
end  if 

print  *,'DO  YOU  WISH  TO  APPLY  A  NON-NEGATIVITY' 
print  *, 'CONSTRAINT?  INPUT  1  FOR  YBS ' 
read  *,non 

print  *, 'INPUT  THB  TYPB  OF  BASIS  FUNCTION  TO  BE  USBD* 

print  * , ' TO  CALCUKARB  THB  UNFOLDED  SPECTRUM* 

print  1-PLANCKIAN  BLACK  BODIES' 

print  *,'  2-CUBIC  SPLINBS' 

read  *,type 

if  (type.eq.l)  then 

print  *,' INPUT  THB  NUMBBR  OF  BASIS  FUNCIONS  TO  BB' 
print  * , ' USED  (MUST  BB  AN  INTEGER  LESS  THAN  15)' 
read  *,nj 

INPUT  THB  INITIAL  GUBSS  AT  THE  PARAMETERS  OF  THE 
BASIS  FUNCTIONS  PROGRAM  ASSUMES  ONE  COBFF .  AND  ONE 
PARAMETER  <T>  PBR  BASIS  FUNCTION 
np»2 . 0*n J 
do  95  j-l.nj 

print  *, 'INPUT  THE  INITIAL  GUESS  FOR  A(J>,J-',J 
read  *,param<  J  > 
continue 
do  100  1-1, nj 

print  *, 'INPUT  THB  INITIAL  GUBSS  FOR  THB  PARAMBTBR 
print  *,'T(J),J-',J, 'TEMP  IS  IN  UNITS  OF  Bo' 
read  *,param<  j+nj  ) 
continue 
end  i  f 

if  (type.eq.2)  then 

print  *, 'INPUT  THE  NUMBBR  OF  KNOTS  TO  BB  USBD  TO  FIT 
print  * , 'THB  CUBIC  SPLINBS  (MUST  BE  AN  INTEGBR' 
print  *, *6<-#KNOTS<-15  )  NOTB:0.0' 
print  *,’IS  COMPUTBD  AS  THB  FIRST  KNOT,  BO( 1  )  IS' 


print  + , 'COMPUTBD  AS  KNOT  2,  AND  128.0*BO<1>  IS* 
print  *, 'COMPUTED  AS  THB  LAST  KNOT.* 
print  *,'INCLUDB  THBSB  KNOTS  IN  THB  COUNT  OP  YOUR' 
print  *, 'KNOTS  AND  BNSURB  THB  KNOTS  ARB  INPUT  PROM' 
print  * , ' THB  MINIMUM  TO  THB  MAXIMUM* 
read  *,nkn 

print  *, 'INPUT  THB  INITIAL  GUBSS  AT  THB  KNOTS' 
print  *, 'KNOT  1-0.0' 
knotl-0.0 
knot2-eo< 1 ) 

print  * , ' KNOT  2-',knot2 
do  101  J-l,nkn-3 

print  *,  'input  knot  *,J  +  2 
read  * ,par am< j  ) 

L01  continue 

knotl5-128 . 0*eo<  1 ) 
print  *, 'KNOT' ,nkn, ,knotl5 
do  102  J-l ,nkn-2 
if  < } .ne. 1 )  then 

print  INPUT  THB  BXPANSION  COBF .  FOR' 
print  *, 'KNOT- ’ ,param(  j-1 > 
else 

print  *, 'INPUT  THB  BXPANSION  COEF .  FOR* 
print  * , ' KNOT- ' , knot2 
end  if 

read  * ,param< nkn-3+ J > 

102  continue 

print  *f 'INPUT  THB  TBMP .  TO  BB  USED  TO  CONSTRUCT  A' 
print  * , *  PLANCKIAN  FIT  BBTWBBN  THB  LAST  2  KNOTS,' 
print  * , ' TEMP  IS  IN  UNITS  OF  Bo' 
read  * ,param( nkn*2-4 > 
np-nkn*2-4 
end  if 

INITIALIZE  THB  H  MATRIX  FOR  THB  MINIMIZATION  TO  THB 
UNITY  MATRIX 
do  110  p*l,np 
do  105  o-l, np 

if  (o.eq.p)  then 
h<  o,p >-1.0 
else 

h<  o,p  >-0.0 
end  if 
105  continue 
110  continue 

INPUT  THB  INITTAL  GUBSS  AT  THB  FUNCTIONAL  VALUE  OF  THB 
LOWBR  BOUND 

print  *, 'INPUT  THB  INITIAL  GUBSS  AT  THB  FUNCTIONAL' 
print  *, 'LOWER  BOUND' 
read  *,1 

BBGIN  THB  APPROXIMATION  OP  THB  TRUB  SPBCTRUM  USING  THB 
UNPOLDBD  SPBCTRUM  FROM  THB  SUBROUTINES  BBFUNC  AND 
CSFUNC 


If  (type.eq.l)  then 

call  bbfunc<  param,su,y,nk,n j ,pi , kc,non  ) 
end  if 

If  <type.eq.2)  then 

call  csfunc<  param,su,nkn,nk,y,non,knotl 
,knot2,knotl5,eo,pi ,kc  ) 
end  1  f 

CALCULATE  C< I  ) 
do  115  1*1, nl 
c<  1  >-0.0 
do  111  k-1 ,nk 

c(  1  )— c(  1  )+su(  k  )*r(  i ,k  )*y 
continue 
continue 
continue 

CALCULATE  CHI  SQUARED 

dummy-chi2 

chl2-0 . 0 

do  120  1-1, nl 

chi2*chi2+(  (  c(  1  >-b<  1  )  )/sigma(  1  >  )**2 
continue 

print  * , ' it- ' , it , *chi2- ' , chi2 
dum3*ab3<  (  dummy-chi2  >/chi2 ) 

If  ( chi2.1e.l .Oe-2 )  then 

print  *, 'minimization  required* , it ,' iterations ' 
print  *, 'chl2- ' ,chl2 
do  122  p-l,np 

print  *, 'parameter- ' ,param(  p ) 
continue 

e-<  eo(  1  )-y/2 .0)+(nk+1.0)*y 
do  130  k-l,nk 
e-e-y 

write  (7,*)  e,  ’  '  ,sa( nk+l-k ) 

continue 

write  (7,*)  eo(l),'  ','0.0' 
if  (type.eq.l)  then 

call  bbfunc<  param,su,y,nk,nj,pi,kc,non  ) 
end  if 

if  <type.eq.2)  then 

call  csfunc( param,su,nkn,nk,y ,non 
'  , knot  1 , knot2, knot  15 , eo, pi , kc ) 

end  if 

e-eo( 1  )-y/2.0 
do  131  k-l,nk 
e-e+y 

write  (7,*)  e,’  ',eu(k) 
continue 
go  to  151 
end  i  f 

if  ( chi2. It .ni .and.dum3 . It . 1 . Oe-2 )  then 

PRINT  ACTUAL  SPBCTRUM  AND  THE  APPROXIMATED  SPBCTRUM 

TO  AN  OUTPUT  PILB 
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print  *, 'minimization  required* , it ,  * iterations 
print  *, *chi2-’ ,chi2 
do  140  p-l,np 

print  * ,  *  parameter- ' ,param(  p ) 
continue 

e*(  eo(  1  )-y/2 .0  >♦<  nk+1  )*y 
do  149  k-l,nk 
e-e-y 

write  (7,*)  e,*  *,sa(nk-k+l> 
continue 

write  ( 7,*  )  eo( 1 >, *  * , *0.0* 
if  (type.eq.l)  then 

call  bbfunc<  param,su, y ,nk,n j , pi,kc,non ) 
end  i  f 

if  <type.eq.2>  then 

call  csfunc<  par am,su,nkn,nk, y,non, 
knotl , knot 2, knot  15 , eo,pi ,kc ) 
end  if 

e*eo<  1  )-y/2 . 0 
do  150  k-l,nk 
e-e+y 

write  (7,*)  e,’  *,su<k) 

continue 
else 

it-it+1 

call  minimizeC  param,y,nk,nj ,pi , kc, non, c,b, sigma 
,np,ni ,h,r , type, nkn, knotl, knot2, knot 15, eo ) 
go  to  118 
end  if 
continue 
close  <  7  ) 
stop 
end 


PLANCKIAN  BLACK  BODY  BASIS  FUNCTION  SUBROUTINE 
subroutine  bbf unc<  param,su,y,nk,nj ,pi ,kc,non  > 
real* 8  eo(  20  ),x,pi,kc,param(  30),su(  1280  ),e,y 
integer  k, J ,nk,n j ,non 
common  eo<  20  ) 
print  *,eo( 1 > 
e°eo( 1  )-( y/2.0  ) 
do  2000  k«l,nk 
e-e+y 
su<  k  )-0.0 
do  1900  3-1, nj 

x-param(  j  >*<  15.0/<  pi*kc*param<  j+n  j  )  >**4  ) 
if  (  e/< kc*param( j+nj > > . gt . 86 . 0 >  then 
su<  k  )-0.0 
else 

su<  k  )-su(  k  >+x*<  e**3.0/<  exp(  e/(  kc*param(  j+nj  ) ) ) 

-1.0  ) ) 


end  1  f 

if  ( non . eq. 1 . and.su< k > . It . 0 . 0 >  then 
su(  k  >-0 . 0 
end  1  f 

1900  continue 
2000  continue 
return 
end 


CUBIC  SPLINE  BASIS  FUNCTION  SUBROUTINE 
subroutine  csfunc<  param,su,nkn,nk,y,non, 

+  knot  1 , knot 2, knot  15 , eo, pi , kc ) 
real+8  param<  30  >,kn( 15  >,coef  (  15  ),nuro/den 
real* 8  planck, cons t ,  pi ,kc , t , x , slope , su< 1280  >, y, e 
real*8  knotl,knot2,knotl5,eo(  20 > 
integer  3,k,non,nk,nkn 
kn( 1 >-knotl 
kn<  2  )-knot2 
do  2010  J*=3,nkn-1 
kn(  3  )-param<  J~2  ) 

2010  continue 

kn<  nkn  )*=knotl5 
do  2020  J-2,nkn-l 

coef <  J  )-param(  J+nkn-4 ) 

2020  continue 

t-param<  2*nkn-4 ) 
e«eo( 1 )-y/2.0 
do  2035  k-l,nk 
e-e+y 
su( k  )-0.0 
do  2030  J-2,nkn 
if  <  j . eq. 2 )  then 

if  (  e .ge .kn(  3  ) . and. e. It .kn( 4 > >  then 

numm( e-kn<  J  +  l  )  )*<  e-kn(  j  +  2 ) )*<  e-kn(  J  +  3 ) > 
den-(  kn(  J  )-kn<  3  +  1  )  )*<  kn(  j  )-kn(  J+2  )  )* 

*  < kn(  j  )-kn( 3  +  3  ) ) 

su(  k )-coef  (  3 )*num/den+su(  k  ) 
end  if 

if  (  e .ge. kn(  2  ) . and . e . It .kn< 3  ) )  then 
slope-1 . 0/<  kn<  2  >-kn( 3  ) ) 

su<  k  )-coef  (  3  )*<  1 .0+<  (  e-kn(  2  >  )*slope  )  >+su(  k  ) 
end  1  f 

if  <  e .  It  .kn<  2  ) . or  . e  .ge .  kn(  4  )  )  then 
su<  k  >-0 . 0+su(  k  ) 
end  i  f 
end  if 

if  (  3 .eq. 3  >  then 

if  < e . ge.kn< 4 ). and.e . It .kn( 5  ) )  then 
if  <nkn.gt.6>  then 

num-<  e-kn<  3  +  1  >  >  +  <  e-kn(  3+2  >  )*<  e-kn<  3  +  3  )  > 


h 

s. 


i 


I 


I 


den-<  kn<  3  )-kn(  3  +  1  >  )*<  kn<  3  >-kn(  3  +  2  )  )* 

*  ( kn( 3  >-kn( 3  +  3  >  > 

su<  k  >-coef <  3  )*num/den+su( k  ) 
else 

su<  k  )-0.0+su<  k  ) 
end  if 
end  If 

if  < e .gt .kn< 3  ) .and.e . It .kn( 4 > >  then 

num-<  e-kn<  3-1  )  >*<  e-kn( 3  +  1 >  )*( e-kn( 3+2  )  ) 
den-<  kn(  3  )-kn<  3-1  >  >*(  kn<  3  >-kn<  3  +  1  ) )+ 

*  <  kn(  3  )-kn<  3  +  2  )  > 

su<  k )-coef (  3 )+num/den+su<  k ) 
end  1  f 

if  (  e.ge . kn(  2 ) . and . e . It . kn( 3 > >  then 
slope-1 . 0/< kn<  3  >-kn<  2  ) ) 
su(  k  )-coef  (  3  >*<  <  e-kn(  2  )  )*slope  )+su(  k  ) 
end  if 

if  < e . It .kn(  2  ) . or . e .ge .kn( 5  )  )  then 
su(  k  )-0 . 0+su<  k  ) 
end  i  f 
end  i  f 

if  < 3 -ge . 4 . and. 3 . It . < nkn-2  ) )  then 

if  < e.ge .kn(  3  +  1 > . and.e . It .kn( 3+2  )  )  then 
if  ( kn< 3+2  ) .ne ,kn< nkn-1  )  )  then 

num-<  e-kn<  J+-1  )  )*<  e-kn<  3+2  )  >*<  e-kn<  3  +  3  >  > 
den-(  kn(  3  )-kn(  3  *1  >  >*<  kn(  3  >-kn(  3  +  2  )  )* 

*  ( kn( 3 >-kn< 3  +  3  ) ) 

su<  k )-coef  (  3 )*num/den+su(  k  ) 
else 

su<  k  )-0 . 0+su(  k  ) 
end  if 
end  if 

if  (  e .ge . kn( 3  ) . and.e . It .kn( 3  +  1 > >  then 
num-(  e-kn<  3-1  )  )*<  e-kn(  3  +  1  >  >*<  e-kn(  3  +  2  ) ) 
den-(  kn(  3  )-kn(  3-1  >  )*(  kn(  3  >-kn(  3  +  1  >  >* 

+  < kn< 3  )-kn< 3+2  ) ) 

su(  k  )-coef(  3  )*num/den+su(  k  ) 
end  i  f 

if  < e .ge .kn( 3-1 > . and . e . It .kn< 3 > >  then 
num-( e-kn<  3-2  )  )  +  ( e-kn<  3-1  )  )  +  ( e-kn( 3  +  1 >  > 
den-<  kn(  3  >-kn<  3-2  )  )*(  kn(  3  >-kn<  3-1  >  >+ 

*  <  kn<  3  )-kn( 3  +  1 > > 

su<  k  )-coef < 3  )*num/den+su( k  ) 
end  i  f 

if  ( e.ge. kn( 3-2  ). and.e. It. kn( 3-1 > >  then 
if  < kn< 3-2 >.ne.kn< 2  )  >  then 

num-<  e-kn<  3-3  >  >*< e-kn<  3-2  )  >+( e-kn( 3-1 > ) 
den-<  kn<  3  >-kn<  3-3  >  >*(  kn<  3  )-kn(  3-2  >  >+ 

*  < kn< 3 >-kn< 3-1 > > 

su(  k  >-coef ( 3  >*num/den+8u( k  ) 
else 

8U<  k  >«0.0+au( k  ) 


r 


end  If 
end  If 

if  < e . It .kn< 3-2 > . or . e . ge. kn(  3  +  2 > )  then 
su(  k  )-0 . 0+su<  k  ) 
end  if 
end  if 

if  ( 3 .eq.< nkn-2 ) )  then 

if  < e .gt .kn< nkn-2 >. and . e . It . kn(  nkn-1 ) >  then 
slope-1 . 0/( kn( nkn-2 )-kn(  nkn-1 ) ) 

su<  k  )“coef  <  3  )*<  1 . 0+<  (  e-kn(  nkn-2  >  >*slope  >)+su(k) 
end  i  f 

if  ( e .ge.kn< nkn-3 >. and . e . It . kn< nkn-2 ) >  then 
num-<  e-kn<  3-2  >  )*<  e-kn<  3-1 >  >*(  e-kn(  3  +  1 >  > 
den-<  kn(  3  )-kn(  3-2  )  >*<  kn(  3  >-kn<  3-1  )  >* 

( kn( 3  >-kn( 3  +  1 >  > 
su<  k  )-coef ( 3  )*num/den+su< k  ) 
end  if 

if  < e .ge.kn< nkn-4  ). and. e  , It .kn( nkn-3  )  )  then 
num-( e-kn(  3-3  >  >*<  e-kn(  3-2 ) >*(  e-kn(  3-1 >  > 
den-(  kn(  3  )-kn(  3-3  ) >*<  kn(  3  >-kn<  3-2  )  )* 

<  kn< 3  )-kn( 3-1  )  > 

su<  k  )-coef < 3  )*num/den+9u( k  ) 
end  if 

if  ( e. It ,kn( nkn-4  ). or .e.ge.kn( nkn-1  ) )  then 
su<  k  )-0.0+su<  k  ) 
end  if 
end  if 

if  < 3 .eq.< nkn-1  ) )  then 

if  < e.ge.kn( nkn-2  ) .and. e. It .kn(  nkn-1 > )  then 
slope-1 . 0/<  kn<  nkn-1 )-kn( nkn-2 ) ) 
su( k  )-coef ( 3  >*( <  e-kn<  nkn-2  )  )*slope  >+su(  k  > 
end  if 

if  ( e .ge . kn( nkn-3  ). and. e . It . kn(  nkn-2 > )  then 
num-<  e-kn( 3-3  )  )*<  e-kn(  3-2 ) >*<  e-kn(  3-1 >  > 
den-(  kn(  3  >-kn(  3-3  >  )*<  kn(  3  >-kn<  3-2  )  >* 

<  kn<  3  >-kn( 3-1 > ) 

su(  k  )-coef(  3  >*num/den+su(  k  ) 
end  if 

if  < e . It .kn( nkn-3  ). or . e .ge . kn(  nkn-1 ) )  then 
su<k)-0.0+su(k) 
end  if 
end  if 

if  (j.eq.nkn)  then 

if  < e.ge.kn( 3-1  ).and.e. le.kn(  3 ) >  then 
x-15 .0/<  pi+kc*t )**4 
planck-x*<  kn<  3-1 >**3.0  >/ 

<  exp(  kn<  3-1 >/( kc*t  >  )-l .  0  > 
const-coef( 3-1 >/planck 

if  ( e/< kc+t > .ge . 69 . 0 )  then 
eu<  k  >-0 . 0+su<  k  > 
else 

su<  k  )-const*x*(  e**3 . 0 >/ 
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*  < exp< e/< kc*t > >-l .0  )+su( k  ) 

end  If 
else 

su<  k  >-0.0+su<  k  ) 
end  If 
end  If 
2030  continue 

If  <non.eq.l.and.su< k). It. 0.0  )  then 
su(  k  >-0.0 
end  if 
2035  continue 
return 
end 
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c  SUBROUTINE  TO  MINIMIZB  CHI  SQUARBD 

subroutine  minimize<  par am, y,nk,n J ,pi ,kc,non,c 
* ,b, sigma, np,nl ,h ,r , t ype,nkn, knot 1 , knot 2, knot 15 ,eo  ) 
real*8  grchi2(  30  ),param<  30  >,b(  20  >,c<  20  ),slgma(  20  ) 
real*8  gradc( 30  ),sprime(  30  >,pprime,t 
real*8  y, pi, kc,r< 20,1280 ),h(  30,30  ),norms ,s( 30  > 
real* 8  v(  30,1  ),gprlme(  30  ),g<  30  ),q(  30,1  ),a(  30,30  ) 
real*8  duml ,sumi ,sum j 

real*8  sumk,su(  1280  ),z<  30,30  ), sump, sumo, pc(  30  >,1 
real*8  kn< 15  ),knotl,knot2,knotl5,eo<  20  > 
integer  nkn,p,np, i ,ni ,o, reset, nk,nj ,non, J ,k,type 
c  CALCULATB  THE  GRADIBNT  OP  CHI  SQUARED 
t-0 

do  2500  p-l,np 
grchi2(  p  )-0.0 

call  gradlentc( param,y,nk,nj , pi ,kc,non,c, 

*  grade, r ,np,ni ,p, type,nkn,knotl , knot 2, knot 15, eo > 
do  2400  i-l,ni 

grchi2(  p  )-grchi2<  p  )+<  (  c(  i  >-b<  i  ) )/<  slgma<  i  >**2  ) ) 

*  *gradc< i  ) 

2400  continue 

grchi2C  p  )-grchi2( p  >*2.0 
2500  continue 
2525  continue 

c  CALCULATB  THB  SB ARCH  DIRECTION 
do  2600  p-l,np 
spr ime( p  )-0.0 
do  2550  o-l, np 

spr  ime(  p  )--h<  p,o  >*grchi2(  o  )*spr  ime<  p  > 

2550  continue 
2600  continue 

c  NORMALIZE  TO  FIND  THB  UNIT  SB ARCH  DIRBCTION 
duml-0.0 
do  2700  p-l,np 

duml-duml+spr lme( p  >**2. 0 
2700  continue 
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norms -sqrt<  duml  > 
do  2800  p-l,np 

s<  p  >-spr ime(  p  >/norms 
2800  continue 

call  llnsearchC np,param, t ,s,y,nk,n j ,pl ,kc,non,ni 
*,r ,b, sigma, h, reset, 1, type, nkn, knot l,knot2,knotl5,eo > 

If  (  reset. eq.l)  then 
reset-0.0 
do  2820  p-l,np 
do  2810  o-l, np 
if  (o.eq.p)  then 
h<  p,o  )-1.0 
else 

h<  p,o  )-0.0 
end  If 
2810  continue 

2820  continue 

go  to  2525 
end  If 

2850  continue 

do  2900  p-l,np 

ppr ime<  p  )-param(  p  )+t*s<  p  > 

2900  continue 

if  <type.eq.2)  then 
do  2903  p-l,nkn-3 
do  2902  J-l,nkn-3 
kn(  j+2)-pprime<  J  ) 

2902  continue 
kn( 1  )-knotl 
kn<  2  )-knot2 
kn( nkn  >-knotl5 

if  < kn< p+2  ). le.kn< p+1 >.or .kn< p+2 >.ge.kn< p+3 > >  then 
t-t/2.0 
go  to  2850 
end  if 

2903  continue 
end  if 

do  2904  p-l,np 

print  *,  'parameter-1  ,param<  p),'s>p)-',s(p>,'t-',t 

2904  continue 

c  CALCULATB  NBW  H  MATRIX  BASBD  ON  NBW  PARAMBTBRS 
do  2905  p-l,np 

v( p, 1 >-ppr ime<  p  >-param<  p  ) 

2905  continue 

c  CALCULATB  A  NBW  C< I >  P0R  NBW  PARAMBTBRS 
if  (type. eq.l)  then 

call  bbfunc<pprime,su,y,nk,nj,pi,kc,non  ) 
end  if 

if  (type.eq.2)  then 

call  csfunc<  ppr ime,su, nkn, nk,y, non, knot 1 
*  ,knot2,knotl5,eo,pi ,kc ) 
end  if 
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do  2908  1-1, ni 
pc<  1  )-0.0 
do  2907  k-l,nk 

pc<  1  >-pc<  i  )+su<  k  >*r<  i  ,k  >*y 

2907  continue 

2908  continue 

do  2915  p-l,np 

call  gradients ppr ime,y,nk, n j ,pi , kc,non, pc 

*  ,gradc,r ,np,ni,p, type,nkn,knotl,knot2,knotl5,eo > 
gpr ime< p  )-0.0 

do  2910  1-1, nl 

gpr  lme<  p  )-gpr  ime(  p  )+< <  pc(  1  >-b<  1 )  >/s igmaC  1  >**2  ) 

*  *gradc< 1 ) 

2910  continue 

gprlme<  p  )-gpr ime( p  >*2 . 0 
2915  continue 

do  2925  p-l,np 

call  gradients  param,y,nk,nj,pi,kc,non,c 

*  , grade, r ,np,ni ,p, type, nkn, knot 1 , knot 2, knot 15, eo  ) 
g(p  )-0.0 

do  2920  1-1, nl 

g<  p  )-g(  p  )+( <  c<  1  >-b<  1  >  >/slgma(  1  )**2  )*gradc(  1  ) 
2920  continue 

g<  p  )«g<  p  )*2.0 
2925  continue 

do  2930  p-l,np 

q<p,l  )-gprime<  p  )-g<  p  ) 

2930  continue 

do  2936  p-l,np 
do  2934  o-l, np 
sumJ-0 .0 
do  2932  1-1,1 

sum j=sumj+v< p, j  )*v( o, j  ) 

2932  continue 

a<  p,o  )-sumJ 
2934  continue 
2936  continue 

do  2942  J-1,1 
do  2940  1-1,1 
sump-0 . 0 
do  2938  p-l,np 

sump-sump+v<  p ,  J  )*q(  p ,  i  ) 

2938  continue 

duml-eump 
2940  continue 
2942  continue 

do  2946  p-l,np 
do  2944  o-l, np 

a<  p,o  )-a<  p,o  )/duml 
2944  continue 
2946  continue 

do  2956  p-l,np 


do  2954  o-l, np 
sum j -0.0 
do  2952  j-l,np 
suml-0.0 
do  2950  1-1,1 
sumk-0.0 
do  2948  k-l,np 

sumk-sumk+q(  k, i >*h<  k,o  > 

2948  continue 

suml-suml+sumk*q<  J,i  ) 

2950  continue 

sumj-sumj+sumi*h<  p, j > 

2952  continue 

z<  p,o  >-sumJ 
2954  continue 
2956  continue 

do  2964  i-1,1 
do  2962  J-1,1 
sump-0.0 
do  2960  p-l,np 
sumo-0 . 0 
do  2958  o-l, np 

sumo-sumo+h<  p,o  )*q<  o, )  ) 

2958  continue 

sump-sump+sumo*q<  p, i > 

2960  continue 

duml-sump 
2962  continue 
2964  continue 

do  2995  p-l,np 
do  2993  o-l, np 

z( p,o  )-z<  p,o  >/duml 

h<  p,o  >-h(  p,o)+a(pro  >-z<  p,o  ) 

2993  continue 
2995  continue 

LBT  PPRIMB  BQUAL  THB  NEW  PARAMBTBRS  AND  RB-BVALUATB  CHI 
SQUARED 

do  3000  p-l,np 

par am( p  )-ppr ime<  p  > 

3000  continue 

LBT  C< I  )-PC< I > 
do  3010  1-1, nl 
c<  i  )-pc(  i  ) 

3010  continue 
return 
end 


SUBROUTINE  TO  CALCULATE  THB  STBP  SIZB  IN  THB  UNIT 
SB ARCH  DIRECTION 

subroutine  llnsearch(  np,param,t,s,y,nk,nj,pl ,kc,non 
*,ni ,r , b, sigma, h, reset , 1 , type,nkn,knotl ,knot2,knotl5,eo > 
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real *8  h( 30,30  ),h2,ht,m,t,to,ho,mo,t2,l,m2,z,duml 
real*8  tl,w,mu,sigma(  20 ) 

real* 8  param( 30  >,s( 30  >,y,pi,kc,r<  20,1280  >,b<  20  > 
real*8  knotl,knot2,knotl5,eo(  20 ) 

Integer  nkn, type, np, reset, nl,nj ,nk,non 
t-0.0 

call  calcfunc<  np,param, t ,s , y ,nk,n ) ,pl ,kc,non 
*,ni,r ,ht, b, sigma, m, type, nkn, knot 1, knot 2, knot 15, eo  > 
ho-ht 
mo-m 

if  (mo.gt.0.0)  then 
reset-1 
return 
end  if 
to-0.0 

if ( ho. le. 1 )  then 
1-ho+C  mo/2.0 ) 
end  if 

if  ( 1.0. It .< -2.0*(ho-l >/mo  ) )  then 
t2-1.0 
else 

t2--2.0*<  ho-1 >/mo 
end  1  f 

3600  continue 
t-t2 

call  calcfunc<  np,par am, t ,s , y ,nk,nj ,pi ,kc,non 
*,ni,r,ht,b, sigma , m , t  ype , nkn , knot 1 , knot 2 , knot 1 5 , eo ) 
h2-ht 
m2-m 

if  (m2.gt.0.0)  then 

z-3 . 0*(  ho-h2  )/<  t2-to  >+mo+m2 

duml-z*2 . 0+mo+m2 

if < duml .eq.0.0  )  then 

tl-to+mo*<  t2-to  >/<  2.0*<  z+mo ) ) 
else 

w-sqrt<  z**2-mo*m2 ) 
mu-<  m2+w-z  >/<  2 . 0*w+m2-mo  > 
tl-t2-mu*<  t2-to ) 
end  if 
end  if 

if < m2. le.0 .0. and.h2 .gt .ho  )  then 
t2-< to+t2 >/2 . 0 
go  to  3600 
end  if 

if  (m2. eq.0.0. and. h2.lt .ho  )  then 
tl-t2 
end  if 

if  ( m2. It . 0 . 0. and.h2 . le.ho )  then 
duml-t2-to 
to-t2 
ho-h2 
mo-m2 
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t2-t2+2.0*<  duml > 
go  to  3600 
end  11 

If  (ro2.eq.0.0.and.ho.eq.h2  >  then 
tl-to 
end  11 
t-tl 
return 
end 


SUBROUTINE  TO  CALCULATE  THB  FUNCTION  TO  MINIMIZB 
H-F< PARAM+T*S> 

subroutine  calc£unc( np,param, t ,s , y , nk,nj ,pi ,  kc 
*,non,nl ,r ,ht , b, sigma ,ra, type, nkn, knot  1 , knot 2, knot 15, eo ) 
real*8  y,tparam<  30  ),param<  30  >,su<  1280  >,b<  20  >,tc(  30  ) 
real*8  t,gradc( 30  ),m,s( 30  >,e,ht ,slgma<  20  > 
real* 8  grchi2<  30  ),pi,kc,r(  20,1280) 
real*8  knotl,knot2,knotl5,eo< 20 > 

Integer  nkn, type,p, i ,  np,ni ,nk,k,non,n j 
do  4000  p-l,np 
tpar am( p  )-param<  p  >+t*s<  p  > 

4000  continue 

CALCULATE  THB  UNFOLDBD  SPBCTRUM 
if  (type.eq.l)  then 

call  bbfunc<  tpar am, su,y,nk,nj ,pi ,kc,non ) 
end  if 

if  <type.eq.2)  then 

call  csfunct  tparam,su,nkn,nk,y,non,knotl 

*  ,knot2,knotl5,eo,pi ,kc ) 
end  if 

CALCULATE  C< I  )  FOR  NBW  PARAMETERS 
do  4200  1-1, ni 
tc( 1  )-0.0 
do  4100  k-l,nk 

tc(  1  >-tc<  i  >+su(  k)*r<l,k>*y 
4100  continue 
4200  continue 

CALCULATE  CHI  SQUARBD  FOR  THB  NBW  PARAMBTBRS 
ht-0.0 

do  4300  1-1, ni 

ht-ht+<  <  tc<  1  )-b(  1  >  >/sigma<  1  )  )**2.0 
4300  continue 

CALCULATE  THB  SLOPB  OF  CHI  SQUARBD 
do  4500  p-l,np 

call  gradients  tpar am, y,nk,n j ,pi ,kc,non 

*  ,tc,gradc,r ,np,nl ,p, type, nkn, knot 1, knot 2, knot 15, eo  ) 
grch!2< p  )-0.0 

do  4400  1-1, nl 

grchl2<  p  )-grchl2<  p  )+<  <  tc<  1  )-b(  1  )  )/slgma<  1  >**2.0  > 

*  *gradc< 1 > 
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4400  continue 

grchl2<  p  >-2.0*grchi2<  p  > 
4500  continue 
m-0.0 

do  4600  p-l,np 

m-m+grchi2<  p  )*s<  p  ) 

4600  continue 
return 
end 


c  SUBROUT IN B  TO  CALCULATE  THB  GRADIENT  OP  C( I  ) 

subroutine  gradients  par am, y ,nk,n J ,pi ,kc,non, c, grade 
*,r ,np,nl,x, type, nkn, knot l,knot2, knot 15, eo > 
real* 6  delta,newpar<  30  >,param<  30  ), f ( 30 >,su< 1280  ) 
real*8  r<  20,1280  >,gradc<  30  ),c<  20  ) 
real* 8  knot 1, Knot 2, knot 15, eo( 20  >,pi ,y,kc 
integer  nkn, i ,x,ni ,p,np,k,nk,n j , non, type 
delta-0.01 
do  5000  p«l,np 
if  ( p.eq.x  )  then 

newparC  p  )-param( p  )+delta*par am< p  > 
else 

newpar<  p  )-param<  p  ) 
end  if 
5000  continue 

if  (type.eq.l)  then 

call  bbfunc(newpar,su,y,nk,nj,pi,kc,non  ) 
end  if 

if  (type.eq.2)  then 

call  csfunc<  newpar ,su, nkn, nk,y, non, knot 1 
*  ,knot2,knotl5,eo,pi ,kc  ) 
end  if 

c  CALCULATE  C< I  )  FOR  THB  NBW  PARAMBTBRS  (  LBT  IT  BQUAL 
c  P(  I  ) ) 

do  5200  1-1, nl 
f<  i  >-0.0 
do  5100  k-l,nk 

f <  1  )-f  (  i  )+su<  k  >*r(  i ,k  )*y 
5100  continue 
5200  continue 

do  5300  1-1, ni 

gradc<  i  >-<  f  <  i  >-c<  i  >  >/<  delta*param<  x  >  > 
if  ( abs< f< i >-c( i > >.lt.le-6 >  then 

print  *, ’WARNING  P( I >-C< I >  LBSS  THAN  1.0B-6’ 
end  if 
5300  continue 
return 
end 
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Appendix  E:  Validation  Results 

This  appendix  is  a  continuation  of  the  results  presented  in  Section  V.  The 
cases  studied  using  one  planckian  basis  function  are  presented  in  Table  XII  and 
the  final  four  cases  studied  using  cubic  spline  basis  functions  with  fixed  knots  are 
presented  in  Table  XIII.  Additionally,  the  linear  plots  of  the  actual  and  unfolded 
spectra  for  the  BC  cases  are  presented.  The  actual  and  unfolded  spectra  for  the 
BP  and  BFI,  BF2,  BF3,  and  BF4  were  identical  so  no  plot  is  presented.  Also,  the 
plots  for  cases  BF5,  BF6,  BF7,  and  BF8  are  presented  even  though  x2  was  greater 
than  20.0. 


TABLE  XII 

Continuation  of  the  Planckian  Basis 
Function  Validation  Test  (a,b) 
Parameters 


Initial 


Unfolded 


a 

T 

a 

T 

2.0 

5.0 

2.0 

5.0 

3.0 

3.0 

2.0 

5.0 

3.0 

4.0 

2.0 

5.0 

1.0 

6.0 

2.0 

5.0 

0.5 

6.0 

2.0 

5.0 

1.0 

2.0 

2.0 

5.0 

1.0 

4.0 

2.0 

5.0 

4.0 

7.0 

2.0 

5.0 

Actual  parameters  equal  2.0  and  5.0 

cr,  =  0.01  and  x2^  0-01  in  all  cases  The  initial  case  is  the  initial  guess  at  the 
unfolded  spectrum 


TABLE  Xin 

Continuation  of  the  Validation  Cases  for  the  Cubic 
Spline  Basis  Functions  With  Fixed  Knots  (a,b,c) 


Case 

Spectra 

Knots 

Coefficients 

T 

r 

BF5 

Actual 

10.0 

25.0 

4.0 

5.0 

3.0 

10.0 

50.0 

80.0 

1.6 

1.0 

Initial 

5.0 

30.0 

2.0 

7.0 

1.0 

2.0 

40.0 

60.0 

3.0 

4,0 

Unfolded 

5.0 

30.0 

3.9 

4.6 

2.5 

14.0 

33.0 

40.0 

60.0 

1.9 

1.5 

BF6 

Actual 

5.0 

10.0 

4.0 

5.0 

3.0 

10.0 

25.0 

50.0 

1.6 

1.0 

Initial 

7.0 

30.0 

2.0 

7.0 

1.0 

20.0 

70.0 

90.0 

3.0 

4.0 

Unfolded 

7.0 

30.0 

4.6 

3.8 

-0.3 

2.0 

2.3e4 

70.0 

90.0 

-17.0 

-81.0 

BF7 

Actual 

5.0 

10.0 

4.0 

5.0 

3.0 

10.0 

25.0 

50.0 

1.6 

1.0 

Initial 

7.0 

30.0 

2.0 

7.0 

1.0 

2.0 

70.0 

90.0 

3.0 

4.0 

Unfolded 

7.0 

30.0 

5.0 

3.5 

1.5 

-12.0 

1.9e3 

70.0 

90.0 

0.55 

-0.17 

BF8 

Actual 

10.0 

25.0 

4.0 

5.0 

3.0 

10.0 

50.0 

80.0 

1.6 

1.0 

Initial(c) 

5.0 

30.0 

2.0 

7.0 

1.0 

20.0 

40.0 

60.0 

3.0 

4.0 

Unfolded 

5.0 

30.0 

3.9 

4.6 

2.5 

14.0 

33.0 

40.0 

60.0 

1.9 

1.5 

a.  cr-  =  0.01 

b.  Convergence  criteria:  x*^  1.0e-2  for  all  cases  or  0.01%  change  in  \2  for  case 
BF5  and  0.1%  change  in  x2  for  cases  BF6,  BF7,  BF8 

c.  Fixed  knots  at  0.0,1.0,128.0  are  implicit  in  definition  of  these  splines. 
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Figure  20.  Linear  Plot  of  the  Actual  and  Unfolded  Spectra 
Versus  Energy  for  Case  BF8 
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Appendix  F:  Test  Case  Results 

This  appendix  is  a  continuation  of  Section  V  and  presents  the  linear  plots  of 


the  cases  presented  in  Section  V.  If  the  linear  plot  of  the  actual  and  unfolded 
spectra  for  a  given  test  case  is  not  in  this  section  and  the  x2  value  for  the  case  in 
Section  V  is  acceptable  (i.e.  less  than  20.0),  then  either  the  actual  and  unfolded 
spectra  match  exactly  or  the  plots  were  presented  in  Section  V. 
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Figure  23.  Linear  Plot  of  the  Actual  and  Unfolded  Spectra 
Versus  Energy  for  Case  NP02 
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Figure  29.  Linear  Plot  of  the  Actual  and  Unfolded  Spectra 
Versus  Energy  for  Case  NP08 
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Figure  34.  Linear  Plot  of  the  Actual  and  Unfolded  Spectra 
Versus  Energy  for  Case  NC02 


r  f  /  /  j  *  /v 


(Tra?7SSgay j  ^SU9W 

Figure  35.  Linear  Plot  of  the  Actual  and  Unfolded  Spectra 
Versus  Energy  for  Case  NC03 
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Appendix  G:  Figure  of  Merit  Study 

In  order  to  determine  the  accuracy  of  the  approximation,  eight  figures  of 


merit  were  considered.  The  goal  of  this  study  was  to  determine  if  any  of  the 
figures  of  merit  could  be  related  to  x2-  By  relating  the  figure  of  merit  to  x2*  & 
method  for  determining  the  accuracy  of  the  unfolding  procedure  may  be  esta¬ 
blished.  The  eight  figures  of  merit  are  as  follows:  the  average  absolute  error, 
AAE;  the  average  absolute  relative  error,  AARE;  the  root  mean  square  average 
absolute  error,  RMS  AAE;  the  root  mean  square  average  relative  error,  R.MS.ARE; 
the  weighted  average  absolute  error,  WAAE;  the  weighted  average  absolute  rela¬ 
tive  error,  WAARE;  the  weighted  root  mean  square  average  absolute  error, 
WRMSAAE;  and  the  weighted  root  mean  square  average  relative  error, 
WRMSARE.  Equations  (22)  thru  (29)  present  the  figures  of  merit  studied  and  the 
approximations  used  by  LCDR.  Kirk  Mathews  to  produce  a  computer  program  to 
calculate  the  figures  of  merit.  The  program  was  validated  using  case  BCl  and  the 
TK  Solver  mathematics  package. 
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Due  to  the  fact  all  the  energy  bin  widths  were  equally  spaced,  the  average  abso¬ 
lute  error  and  average  relative  error  can  be  calculated  by  dividing  Equations  (22) 
and  (23)  by  nk,  the  number  of  energy  bins  used  to  evaluate  the  integrals. 
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(25) 

(26) 

(27) 


Finally  the  weighted  figures  of  merit  were  calculated  by  weighting  the  absolute 
error  and  relative  error  with  a  factor  of  l/E. 


Thus 


UX4E 


-/ 


|5„(£)-5.(£)| 


£ 


f  1 5,  (£*)-■?«(£*)  I 


ir,t££ 


=  / 


|5.(£)-5a(£)| 


5U(£)+Sa(£) 


^)<i£ 


E 

* 


K  (£*)-5a(£*)  | 


5.(£*)+5a(£*) 


2S 


(29) 


The  weighted  figures  of  merit  can  then  be  calculated  using  Equations  (21)  thru 
(27)  by  replacing  the  absolute  error  and  the  relative  error  with  the  weighted  abso¬ 
lute  error  and  the  weighted  relative  error  respectively. 

The  results  of  this  study  are  presented  in  Table  XTV. 
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TABLB  XIV 


Comparison  Of  Figures  Of  Merit 
For  Cases  Studied  (a) 


RMS 

RMS 

WRMS 

WRMS 

Cs 

X3 

AAB 

ARB 

AB 

RB 

WAAB 

WAARB 

AB 

RB 

NPO 

0.010 

4e-4 

0.012 

7e-4 

0.017 

5e-4 

0.010 

9e-4 

0.010 

NCO 

0.020 

0.017 

0.047 

0.027 

0.066 

0.028 

0.020 

0.042 

0.035 

BC4 

0.026 

0.001 

0.001 

0.003 

0.002 

0.005 

0.002 

0.008 

0.002 

BC2 

5.4 

0.014 

0.013 

0.029 

0.017 

0.037 

0.014 

0.053 

0.018 

NC04 

5.5 

0.06 

0.11 

0.12 

0.13 

0.23 

0.091 

0.34 

0.11 

NP04 

6.0 

0.005 

0.095 

0.008 

0.11 

0.009 

0.086 

0.012 

0.096 

NC07 

7.2 

0.10 

0.28 

0.14 

0.38 

0.11 

0.11 

0.15 

0.20 

NCO  3 

8.6 

0.13 

0.42 

0.19 

0.59 

0.28 

0.19 

0.33 

0.31 

NC06 

8.8 

0.064 

0.051 

0.13 

0.073 

0.23 

0.078 

0.32 

0.10 

NC05 

9.4 

0.099 

0.44 

0.22 

0.67 

0.27 

0.20 

0.39 

0.36 

NP07 

9.7 

0.003 

0.093 

0.004 

0.12 

0.003 

0.036 

0.004 

0.063 

NC01 

11.0 

0.14 

0.35 

0.22 

0.46 

0.35 

0.19 

0.44 

0.26 

NP06 

12.0 

0.004 

0.048 

0.007 

0.051 

0.005 

0.030 

0.008 

0.037 

NPOS 

13.0 

0.006 

0.17 

0.012 

0.23 

0.010 

0.10 

0.015 

0.14 

NC09 

13.0 

0.12 

0.15 

0.18 

0.19 

0.23 

0.11 

0.30 

0.14 

NP09 

15.0 

0.004 

0.18 

0.006 

0.22 

0.007 

0.097 

0.009 

0.13 

NCOS 

15.0 

0.080 

0.065 

0.16 

0.081 

0.29 

0.11 

0.44 

0.14 

NP03 

16.0 

0.002 

0.036 

0.004 

0.38 

0.004 

0.037 

0.005 

0.040 

PC  3 

16.0 

0.005 

0.071 

0.008 

0.22 

0.009 

0.46 

0.012 

0.86 

NPOS 

17.0 

0.005 

0.086 

0.011 

0.10 

0.012 

0.092 

0.017 

0.10 

NC02 

17.0 

0.18 

0.19 

0.31 

0.22 

0.51 

0.19 

0.65 

0.22 

CP3 

17.0 

0.14 

0.23 

0.28 

0.30 

0.57 

0.22 

0.84 

0.30 

CP4 

17.0 

0.14 

0.27 

0.28 

0.45 

0.56 

0.23 

0.84 

0.33 

NC10 

18.0 

0.038 

0.067 

0.073 

0.082 

0.11 

0.049 

0.16 

0.063 

NPO  1 

18.0 

0.004 

0.28 

0.005 

0.42 

0.003 

0.092 

0.004 

0.21 

NP10 

18.0 

0.003 

0.042 

0.005 

0.049 

0.004 

0.027 

0.006 

0.032 

BC3 

19.0 

0.071 

0.13 

0.089 

0.22 

0.094 

0.050 

0.12 

0.11 

CP2 

20.0 

0.12 

0.086 

0.30 

0.12 

0.59 

0.20 

0.90 

0.29 

NP02 

29.0 

0.003 

0.15 

0.004 

0.20 

0.003 

0.054 

0.004 

0.10 

BP  5 

33.0 

0.081 

0.20 

0.93 

0.35 

0.078 

0.064 

0.091 

0.17 

BF8 

33.0 

0.081 

0.20 

0.93 

0.35 

0.078 

0.064 

0.091 

0.17 

PC  2 

38.0 

0.003 

0.058 

0 . 003 

0.14 

0.003 

0.22 

0.004 

0.57 

PCI 

71.0 

0.020 

0.29 

0.029 

0.43 

0.028 

0.77 

0.039 

1 .11 

BC1 

230.0 

0.094 

0.17 

0.12 

0.22 

0.10 

0.096 

0.14 

0.16 

CPI 

380.0 

0.60 

1.7 

0.85 

1.8 

1.3 

1.0 

1.6 

1.3 

BP  7 

1,900 

0.12 

0.79 

0.23 

1.2 

0.37 

0.29 

0.51 

0.61 

BP  6 

23,000 

0.42 

1.3 

0.59 

1.6 

0.50 

0.58 

0.61 

0.94 

a.  All  figures  of  merit  for  case  BP1,  BP2,  BP3,  BP1,  BP2, 
BF3,  and  BF4  were  less  than  0.01 
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Abstract 

The  purpose  of  this  study  was  to  develop  the 
methodology  for  and  to  implement  a  computer  program  to 
approximate  a  solution  to  a  system  of  Fredholm  integral 
equations.  The  system  of  equations  used  in  this  study  is 
representative  of  the  equations  formed  during  the  detection 
of  pulsed  radiation  using  a  series  of  detectors  with 
asymmetric  response  functions.  Though  general  in  nature  and 
applicable  to  all  systems  of  fredholm  integral  equations, 
the  equations  studied  are  of  importance  to  the  Defense 
nuclear  Agency  with  regard  to  the  measurement  of  radiation 
spectra  during  underground  nuclear  effects  simulation 
testing. 

The  deconvolution  technique  consisted  of  representing 
the  unfolded  spectrum  as  a  weighted  sum  of  basis  functions. 
This  unfolded  spectrum,  the  actual  spectrum,  and  a  predicted 
spectrum  were  then  used  to  form  a  X~  test  statistic.  By 
adjusting  the  parameters  in  the  basis  functions  and  their 
weights,  X*-’  was  minimized  and  the  unfolded  spectrum  was 
corrected  to  approximate  the  actual  spectrum. 

The  methodology  for  this  deconvolution  technique  was 
then  converted  into  a  general  computer  program.  The 
validation  cases  conducted  on  the  two  types  of  spectra 
confirmed  the  reliability  of  the  methodology  and  the 
computer  program.  Additionally,  an  initial  study  with 
simulated  measurement  error  added  to  the 

mea3ur ed-to-predicted  ratios  showed  that  the  actual  spectrum 
could  not  be  returned  exactly.  The  second  study 
approximated  the  actual  spectrum  with  an  unfolded  spectrum 
using  a  second  3et  of  basis  functions.  An  acceptable 
approximation  was  conducted;  however,  certain  artifacts  were 
discovered  in  the  unfolded  spectrum.  The  validation  cases 
and  preliminary  test  cases  conducted  prove  that  the  computer 
program  based  on  the  methodology  presented  in  this  study  is 
a  viable  means  of  approximating  an  actual  radiation 
spectrum . 
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